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ABSTRACT 


The  plane  wave  reflection  coefficient  is  an  important  geometry  independent  means  of  speci¬ 
fying  the  acoustic  response  of  a  horizontally  stratified  ocean  bottom.  It  is  an  integral  step  in  the 
inversion  of  acoustic  field  measurements  to  obtain  parameters  of  the  bottom  and  it  is  used  to 
characterize  an  environment  for  purposes  of  acoustic  imaging.  This  thesis  studies  both  the  gen¬ 
eration  of  synthetic  pressure  fields  through  the  plane  wave  reflection  coefficient  and  the  inversion 
of  measured  pressure  fields  to  estimate  the  plane  wave  reflection  coefficient.  These  are  related 
through  the  Sommerfeid  integral  which  is  in  the  form  of  a  Hankel  transform.  The  Hankel 
transform  is  extensively  studied  in  this  thesis  and  both  theoretical  properties  and  numerical 
implementations  are  considered.  These  results  have  broad  applications.  When  we  apply  them  to 
the  generation  of  synthetic  data,  we  obtain  hybrid  numerical-analytical  algorithms  whidi  pro¬ 
vide  extremely  accurate  synthetic  fields  without  sacrifising  computational  speed.  These  algorithms 
can  accurately  incorporate  the  effects  of  trapped  modes  guided  by  slow  speed  layers  in  the  bot¬ 
tom.  We  also  apply  these  tools  to  study  the  inversion  of  measured  pressure  field  data  for  the 
plane  wave  reflection  coefficient.  We  address  practical  issues  associated  with  the  inversion  pro¬ 
cedure  including  removal  of  the  source  field,  sampling,  field  measurements  over  a  finite  range, 
and  uncontrolled  variations  in  source-height.  A  phase  unwrapping  and  associated  interpolation 
scheme  is  developed  to  handle  improperly  spaced  data. 

A  preliminary  inversion  of  real  pressure  field  data  is  performed.  In  parallel,  an  inversion 
of  a  synthetically  generated  field  for  smilar  bottom  parameters  is  also  performed  and  the  results 
of  processing  the  real  and  synthetic  data  are  compared.  The  estimate  for  the  depth  dependent 
Green’s  function  obtained  from  the  real  data  shares  many  features  with  the  depth  dependent 
Green’s  function  **imatrei  from  the  synthetic  data,  suggesting  that  the  total  inversion  to  obtain 
the  plane  wave  reflection  coefficient  will  soon  be  possible.  Errors  in  the  present  estimate  of  the 
plane  wave  reflection  coefficient  are  associated  with  uncontrolled  source-height  variations  during 
the  acquisition  of  data. 
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CHAPTER  I: 

INTRODUCTION 

l.l)  Overview 

The  plane  wave  reflection  coefficient  is  an  important  geometry  independent  means  of  speci¬ 
fying  the  acoustic  response  of  a  horizontally  stratified  ocean  bottom.  It  is  used  both  to  calculate 
fields  in  propagations  models  and  as  an  input  to  a  variety  of  inverse  techniques  which  seek  to 
determine  bottom  parameters  [1,2,3  ].  In  this  thesis  we  will  study  both  the  generation  of  syn¬ 
thetic  pressure  fields  through  the  plane  wave  reflection  coefficient  and  the  inversion  of  measured 
pressure  fields  to  estimate  the  plane  wave  reflection  coefficient.  We  will  consider  only  the  fields 
associated  with  a  CW  point  source  in  the  deep  ocean  over  a  horizontally  stratified  bottom  and 
will  not  allow  the  bottom  to  support  shear  waves.  The  results  of  this  thesis,  however,  are  appli¬ 
cable  to  a  wide  dasi  of  wave  problems  and  can  be  generalized  to  permit  the  source  to  be  within 
any  isovelocity  layer  and  the  introduction  of  shear.  Further,  it  is  our  hope  that  the  tools 
developed  in  the  course  of  this  work  will  find  applications  in  many  areas. 

The  foundation  for  our  studies  of  the  forward  and  inverse  problems  is  the  Hankel 
transform,  [4,  S,  6  ]  which  arises  in  these  contexts  because  the  Sommerfeld  integral  which  relates 
the  plane  wave  reflection  coefficient  to  the  reflected  field  is  in  that  form.  We  will  derive  general 
properties  of  the  Hankel  transform  to  guide  the  work  in  these  areas.  We  will  also  study  and 
develop  numerical  implementations  to  permit  computer  processing.  A  fast,  accurate  numerical 
Hankel  transform  algorithm  is  developed  and  illustrated. 

The  Hankel  transform  results  allow  us  to  isolate  significant  sources  of  degradation  in 
numerically  generated  synthetic  fields.  To  remove  these  sources  we  develop  a  hybrid  analytical- 
numerical  procedure  which  is  significantly  more  accurate  without  sacrificing  computational 
speed.  This  hybrid  algorithm  performs  some  of  the  calculations  analytically  while  keeping  the 
numerical  computation  in  the  form  of  a  Hankel  transform. 

Through  another  hybrid  technique  we  incorporate  the  effects  of  trapped  modes  that  may  be 
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guided  by  low  speed  layers  in  the  bottom.  The  results  are  accurate  both  in  the  near  and  far 
fields,  in  contrast  with  modal  methods.  The  method  developed  for  handling  these  modes  can 
also  be  used  to  control  other  complications  associated  with  poles  in  the  plane  wave  reflection 
coefficient  such  as  those  that  would  be  introduced  by  allowing  for  shear  wave  propagation. 

In  Chapter  V  we  begin  to  study  the  inversion  of  pressure  field  data  to  obtain  an  estimate 
for  the  plane  wave  reflection  coefficient.  We  draw  upon  our  previous  results  to  consider  a 
recently  proposed  method  for  performing  this  inversion  by  Frisk,  Oppenheim  and  Martinez  [7  ]. 
Frisk,  Oppenheim  and  Martinez  have  proposed  that  the  Sommerfeld  integral  be  inverted 
directly,  without  recourse  to  the  specular  approximation  used  by  previous  methods.  [7  ]  With 
such  a  direct  inversion,  estimates  would  no  longer  be  confined  to  real  angles  of  incidence  and 
regions  where  the  specular  approximation  is  valid.  [8, 9  ]  Such  a  direct  inversion  has  been  made 
possible  recently  by  coherent  measurements  of  the  reflected  pressure  field  resulting  from  a  point 
source  over  a  horizontally  stratified  bottom.  [10  ]  In  this  chapter  we  study  several  practical  issues 
associated  with  this  proposed  direct  inversion.  We  consider  directly  the  issues  of  source  field  sub¬ 
traction,  sampling,  windowing  of  the  pressure  field,  and  uncontrolled  variations  in  source  height. 
The  issue  of  source-field  subtraction  arises  because  the  plane  wave  reflection  coefficient  is 
directly  related  to  the  reflected  field  and  not  the  total  field,  which  is  measured.  Under  the  issue 
of  sampling  we  study  both  the  sampling  rate  required  to  obtain  a  valid  inversion  and  the  effects 
of  improperly  spaced  data.  To  handle  improperly  spaced  data,  an  interpolation  procedure  is 
developed  which  is  based  on  a  new  phase  unwrapping  procedure.  In  the  discussion  of  window¬ 
ing  we  determine  the  range  to  which  field  measurements  must  be  taken  in  order  to  obtain  a  valid 
inversion.  In  the  section  on  source-height  variation  we  characterize  the  degradation  that  results 
from  variation  in  the  source-height  and  study  canonical  variations. 

Having  conadered  many  of  the  issues  affecting  the  direct  inversion  of  pressure  field  data  to 
obtain  the  plane  wave  reflection  coefficient  we  now  actually  perform  a  preliminary  inversion  of 
real  data.  In  parallel  we  invert  synthetic  data  that  we  have  generated  using  bottom  parameters 
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comparable  to  those  we  believe  associated  with  the  real  data.  We  draw  upon  the  previous 
developments  to  interpret  the  results. 

In  the  remaining  portion  of  this  chapter  we  briefly  develop  the  acoustic  framework  upon 
which  this  thesis  rests.  We  also  describe  the  experimental  paradigm  by  which  the  real  data  was 
gathered. 


1.2)  Plane  Waves  and  a  Horizontally  Stratified  Environment 

A  horizontally  stratified  environment  is  one  for  which  the  material  parameters  vary  only 
vertically.  Such  a  simple  environment  makes  possible  an  in  depth  study  without  the  complications 
that  a  more  varied  environment  would  introduce.  The  insights  gained  from  studying  this  simple 
environment  can  provide  an  understanding  of  more  complicated  envonments.  Also,  for  many 
conditions,  such  as  are  found  in  the  region  of  an  abyssal  plane  in  the  deep  ocean,  the  model  is 
itself  sufficient. 

Figure  1.2.1  shows  an  isovelocity  water  layer  over  a  horizontally  stratified  bottom.  Within 
the  water  a  single  plane  wave  has  the  form  t1 

el(k,x+k,y+k,i)e-iul 

For  this  wave  to  be  a  solution  to  the  wave  equation  within  the  water  (which  has  sound  speed  c  ): 

v2~jr$;  *(*’0  =  °  (2) 

the  wave  numbers  (kx ,  ky ,  and  kt)  must  satisfy: 

*x2+*;+*,2  -4  =  o  (3) 

c 

We  define  the  vector  k  »  kxix  +kyiy  +ktit  and  the  scalar,  k  *  — .  k  will  point  in  the  direc¬ 
tion  that  the  plane  wave  propagates.  In  terms  of  this  vector,  the  requirement  (3)  can  be  written 


1)  Throughout  tbli  theas  we  will  suppress  the  e 
the  field  expressions. 


time  dependence  because  it  separates  out  from  sD 


INCIDENT 

PLANE 

WAVE 
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If  (4)  is  satisfied  then  (1)  is  formally  a  solution  to  the  wave  equation  even  if  k  is  complex.  When 
one  or  more  of  the  components  of  k  are  imaginary  the  plane  wave  is  called  evanescent  and  will 
vary  exponentially  in  the  direction  of  the  imaginary  component(s). 


Considerations  of  symmetry  guarantee  that  when  a  plane  wave  strikes  a  horizontally  strati¬ 
fied  bottom  the  resulting  wave  will  also  be  planar.  For  the  purposes  of  field  calculations,  plane 
waves  are  eigenfunctions  of  horizontally  stratified  systems.  An  incident  plane  wave  given  by  : 


Ptel^k,x  +k>}  +i,,) 

will  generate  the  reflected  wave 


(5) 


PKe‘(k,x+l,y  k,,)  (6) 

The  change  in  sign  indicates  that  the  reflected  wave  is  returning  in  the  z  direction.  The  ampli¬ 
tude  change  defines  the  plane  wave  reflection  coefficient,  whidi  may  be  a  function  of  k .  Since 
for  fixed  u>  only  two  of  the  three  components  of  k  can  be  specified  independently  we  write  the 
plane  wave  reflection  coefficient  explicitly  as  a  function  of  only  two: 


rc(**.M»7r  (?) 

Our  plane  wave  reflection  coefficient  is  defined  for  a  single  frequency  only.  This  implies 
that  the  incident  signal  has  been  present  for  all  time.  The  returning  wave,  PKel^k‘x*k,'>  is 
influenced  by  the  bottom  at  all  depths.  This  is  in  contrast  to  the  occasional  usage  of  the  term  for 
which  only  the  surface  contribution  to  a  pulsed  input  is  considered. 


1.3)  The  Weyl  Integral 

Because  the  wave  propagation  we  consider  is  linear  we  can  determine  the  response  to  a 
more  complicated  incident  field  by  considering  that  field  as  a  superposition  of  plane  waves.  [9  ] 

To  calculate  the  reflected  responses  to  a  point  source  shown  in  Figure  1.3. 1  we  first  express 
the  field  of  the  point  source  at  z  =  z0  as  a  superposition  of  plane  waves.  We  write:1 
1)  We  will  uie  to  denote  the  wave  number  in  the  water. 
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Ftgare  1.3.1  Point  source  geometry 
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r,{x,y,*o)  =  fA(kx,k,)elll‘J+k’y)dkxdk,  (i) 

y/x2+y2  Z-n  » 

Equation  (1)  is  a  two  dimensional  Fourier  transform  and  can  be  inverted  to  determine 


*«■*>  -  ■  vrro  (2) 

_  ,  .  .  ...  iVki-k,1- kMt) 

Each  plane  wave  m  the  superposition  (1)  propagates  m  the  z  direction  as  t  or 

e  depending  on  its  direction.  Because  the  incident  waves  are  those  which  carry 

energy  from  the  source  to  the  bottom  we  know  that  all  the  waves  have  Re  (kt)  >  0. 1  This  allows 
us  to  write  the  field  at  z  >  zg  as 


*>  zq  P,(x,y,z)  =  j-]  /  rj-™1—  1elVk*~'^')el(k'x+t>y)dkIdky  (3) 
2-rr  Vki  -kj-kj 

When  these  waves  strike  the  bottom  at  z  =0  each  is  turned  around  (+i,  -  ~kx)  and  is 
scaled  by  the  plane  wave  reflection  coefficient.  The  reflected  field  at  z  =0  is  given  by: 

»  /.  m  Iff  »rc(*x>*r)  .  iVki  -i’-k.H I„) .  i(k,x  +*,y ) 

-  uUvtFJ-Z  ’ 

The  reflected  field  at  any  height  z  >0  is  determined  by  propagating  this  reflected  field  up 
just  as  we  propagated  the  incident  field  down. 

*+*o>0  PK(x,y,z)  =  ~ff  £&&JtLe,Vki  +^el{k‘J  +k>y)dkxdk,  (5) 


Equation  (S)  is  often  called  the  Weyl  integral.  [13  ] 


1.4)  The  Sommerfeid  Integral 


1)  The  oompVra  Poyntiaa__MCtDl-  [11  }  modeled  with  ■  plane  wen  ii  given  by 

I  Vt)  — t J— jLl 

S(x  ct  j)  m  ^PV’  »  - — -2-  \P  I*.  [12  1  The  power  of  the  wive  Qowi  in  the  direction  of  the 

real  pert  of  die  Pointing  vector  A  wire  with  poritive  red  pert  of  cerriei  energy  in  the  po- 

dtiver  direction. 
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The  field  associated  with  a  point  source  over  a  horizontally  stratified  bottom  is  symmetric 
about  the  z  axis  of  Figure  1.3.1.  We  can  exploit  this  symmetry  to  reduce  the  two  dimensional 
Weyl  integral  (1.3.5)  to  a  one  dimensional  integral. 

The  symmetry  of  the  problem  guarantees  that  all  the  field  variables  in  space  show  a  cylindr¬ 
ical  symmetry.  Because  the  two  dimensional  Fourier  transform  of  a  circularly  symmetric  function 
will  also  be  circularly  symmetric,  die  Fourier  domain  will  also  display  a  cylindrical  symmetry. 
We  define 

r 2  ■  x2+y2 

With  these  definitions  we  write  (1.3. 5)  in  cylindrical  coordinates  as: 


_  2lt  oo  V*o -k2 

With  TiVkt+k?)  m  rc(tx  ,1:^).  Performing  the  0  integration  and  using  (14  ] 


w  -  ^~felxc°*d9 

Equation  (2)  becomes 


(2) 


(3) 


P*(r,z)  =  f—7  [  -—r(tf)e'V*^(l^,)/o(trr)trdStf 
o  \ki-k; 


This  is  the  Sommerfdd  integral.  [13  ] 

The  Sommerfdd  integral  has  the  form 


P*(r  .*)  =  fG(kr,xj0)Jo(rkr)krdkr  (5) 

0 

where 

G  (kr  ,x  ,xq)  -  )e,V*^'*+,)  (6) 

We  will  refer  to  G(kr,z ,z0)  as  the  depth  dependent  Green’s  function  or  the  Green's  function. 
The  integral  transform  (5)  is  the  Hankel  transform.  [4, 5, 6  ] 
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Equation  (4)  represents  the  reflected  field  as  a  superposition  of  cylindrical  waves  of  the 
form  J(>(krr)el  each  of  which  can  be  considered  to  be  a  superposition  of  plane 

waves  all  striking  the  surface  with  the  same  horizontal  wave  number,  kr.  Because  F(kr)  is 
related  so  directly  to  the  plane  wave  reflection  coefficient,  V(\/kJ+kJ)  =  Tc(kx  ,k}),  we  will 
refer  to  it  as  the  plane  wave  reflection  coefficient. 


1.5)  Obtaining  the  Reflection  Coefficient 

Presently  most  techniques  for  determining  the  plane  wave  reflection  coefficient  as  a  func¬ 
tion  of  horizontal  wave  number,  r(Jt, ),  from  the  reflected  pressure  field,  Pg  (r  ),  do  not  concen¬ 
trate  on  inverting  Equation  (1.4.4).  Instead  they  consider  the  stationary  phase  approximation  to 
that  Equation  given  by:1  [9  ] 


<■.(')  -  — - 


(1) 


where  &0  i*  defined  to  be  the  water  wave  number,  — ,  and  R2  *  r2+(z0+z)2 

c 


Equation  (1.4.4)  is  inverted  to  provide 


“0* 


P*(r) 


(2) 


fh  fk 

If  Equation  (2)  is  used  to  estimate  T(— -)  then  |r(— -)  |  can  be  estimated  from  the  mag- 

K  A 

nitude  of  Pr  (r)  alone  through: 


|IX^)I-  1*1  \p(r)\ 

This  fact  together  with  the  simplicity  of  Equation  (2)  account  for  its  widespread  use. 


(4) 


1)  Mora  frequently  the  plane  wave  reflection  coefficient  u  a  function  of  angle  is  related  through  the  •pett¬ 
ier  approximation  to  reflected  pr  enure  field.  We  preeent  it  as  ■  function  of  horixonml  wave  number  to  be 
consistent  with  the  ram  of  the  test  If  we  denote  1»(8)  as  the  plena  wave  reflection  ooeflSdeet  is  a  Amo¬ 
tion  of  ingle,  the  two  formulations  an  related  through  r,(B)  =■=  r(Jto**n(0)).  iojj"  1*  *8*  borijootal 
wive  number  correepondlag  to  the  specular  angle. 
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Unfortunately  the  stationary  phase  approximation  upon  which  (4)  is  based  is  appropriate  only 
for  distances,  R ,  large  compared  to  a  wavelength  and  only  for  specular  angles  less  than  criti- 
cal.[8,9]  It  completely  ignores  near  field  effects  associated  with  T(kr)  for  kr  >  k0.  For  applica¬ 
tions  that  consider  near  field  effects  or  the  character  of  the  pressure  field  close  to  or  greater  than 
the  critical  angle,  a  more  exact  inversion  of  Equation  (1-4.4)  is  required. 

For  such  applications  Frisk,  Oppenheim,  and  Martinez  [7]  have  proposed  that  both  the 
magnitude  and  phase  of  the  reflected  pressure  field  be  measured  and  that  the  Sommerfeld 
integral  of  Equation  (1.4.4)  be  inverted  directly,  without  recourse  to  the  stationary  phase 
approximation.  The  Sommerfeld  integral  is  in  the  form  of  a  Hankel  transform.  Since  the 
Hankel  transform  is  its  own  inverse  (5],  Equation  (1.4.4)  can  be  inverted  and  solved  for  the 
plane  wave  reflection  coefficient  in  terms  of  Pp(r),  the  reflected  response  to  a  point  source. 
This  gives: 

IX*,)  =  -iVki-k}e-lV^  U+',l/P,(ry0(MWr  (5) 

0 


1.6)  Experimental  Model 

In  this  thesis  we  will  perform  a  preliminary  inversion  of  measured  pressure  field  data 
according  to  Equation  (1.5.5),  as  suggested  by  Frisk,  Oppenheim,  and  Martinez.  [7]  The  data 
we  analyze  was  taken  by  G.  Frisk,  J.  Doutt,  and  E.  Hays  in  1961 .  In  this  section  we  present  the 
details  of  the  experimental  configuration  they  used.  A  similar  experimental  configuration  has 
been  described  in  the  literature.  [10] 

Figure  1.6.1  shows  the  experimental  configuration  used  by  Frisk,  Doutt,  and  Hays.  As 
shown,  two  receivers  were  moored  1.17  and  54.55  meters  from  the  bottom  of  the  ocean  on  an 
abyssal  plain  under  1900  meters  of  ocean.  The  source  was  attached  by  cable  to  a  lhip  on  the  sur¬ 
face,  and  drifted  slowly  away  at  a  height  off  the  bottom  of  approximately  135  meters.  Every  12 
seconds  die  source  emitted  a  4  second  220  Hertz  tone  which  die  receivers  recorded  after  quadra- 
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Fignre  1.6.1  Experimental  configuration  used  by  Frisk,  Doutt  and  Hays  to  obtain  real  data 
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ture  demodulation  and  low  pass  filtering.  In  this  way  the  complex  amplitude  of  the  field  as  a 
function  of  range  was  recorded.  The  source  strength  was  177  dB  re  1  p.Pa  at  1  meter. 

Recording  by  the  receivers  was  initiated  every  12  seconds  upon  recognition  of  an  11  kHz 
trigger  pulse  sent  from  a  pinger  mounted  on  the  source,  and  continued  for  6  seconds.  During  this 
time  the  output  of  the  receivers  was  quadrature  demodulated,  low  pass  filtered  to  2  Hz,  digitized 
by  a  12  bit  A/D  converter  at  a  5  Hz  rate  and  recorded  on  cassette  tape.  A  schematic  of  the  prel¬ 
iminary  data  processing  in  the  receiver  system  is  shown  in  Figure  1.6.2. 

The  ship  drifted  at  a  speed  of  about  1/2  kn  allowing  one  sample  of  the  field  every  half 
wavelength.  The  clocks  in  the  source  and  receiver  were  synchronized  and  had  a  stability  of  about 
one  part  in  109  per  day.  The  11  kHz  emission  times  at  the  source  and  the  arrival  time  at  the 
receivers  were  used  to  determine  the  slant  range  between  the  source  and  the  receivers.  As  part  of 
die  processing  it  was  necessary  to  estimate  the  source  height  and  convert  from  slant  ranged  to 
horizontal  range. 

Frisk,  Doutt  and  Hays  determined  that  the  signal  was  in  a  steady  state  condition  by  the  4th 
date  sample. 
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1.7)  Summary 

In  this  thesis  we  consider  the  generation  of  synthetic  pressure  fields  through  the  evaluation 
of  the  Sommerfeld  integral.  This  integral  is  in  the  form  of  the  Hank  el  transform  of  the  depth- 
dependent  Green’s  function.  We  also  consider  the  inversion  of  a  measured  pressure  field  to  esti¬ 
mate  the  depth-dependent  Green's  function  and  from  that  the  plane  wave  reflection  coefficient. 
The  foundation  of  both  these  procedures  is  the  Hankel  transform.  In  the  next  chapter  we  will 
both  catalogue  and  develop  the  properties  of  the  Hankel  transform  that  will  provide  the  founda¬ 
tion  for  the  work  of  this  thesis. 


i 


9 


t 


\ 

1 


-  23  - 


References 

1.  D.C.  Stickler  and  P.A.  Drift,  “Inverse  Problem  for  a  Stratefied  Ocean  and  Bottom,"  J. 
Acoust.  Soc.  Am.  70,  pp. 1723-1727  (1981). 

2.  D.C.  Stickler,  Personal  communication. 

3.  H.E.  Moses  and  C.M.  dcRidder,  Properties  of  Dielectrics  from  Reflection  Coefficients  in 
One  Dimension,  Lincoln  Laboratory  Technical  Report  no.  322,  Lincoln,  Massachusetts 
(July  11  1963). 

4.  G.  N.  Watson,  Theory  of  Bessel  Functions  2nd  Ed.,  Cambridge  at  the  University  Press,  New 
York  (1966). 

5.  R.  Bracewell,  The  Fourier  Transform  and  Its  Applications,  McGraw-Hill,  New  York  (1965). 

6.  Papoulis,  Systems  and  Transforms  with  Applications  in  Optics,  McGraw-Hill,  New  York 
(1968). 

7.  G.V.  Frisk,  A.V.  Oppenheim,  and  D.  Martinez,  “A  Technique  for  Measuring  the  Plane- 
Wave  Reflection  Coefficient  of  the  Ocean  Bottom,"  J.  Acoustical  Sot.  Amer.  68  (2), 
pp. 602-612  (1980). 

8.  D.C.  Stickler,  “Negative  Bottom  Loss,  Critical-Angle  Shift,  and  the  Interpretation  of  the 
Bottom  Reflection  Coefficient,”  J.  Acoust.  Soc.  Am.  61,  pp.707-710  (March  1977). 

9.  Brekhovsldkh,  Waves  in  Layered  Media,  Academic  Press,  New  York  (1960). 

10.  G.V.  Frisk,  J.  Doutt,  and  E.  Hays,  “Bottom  Interaction  of  Low-Frequency  Acoustic  Sig¬ 
nals  at  Small  Grazing  Angles  in  the  Deep  Ocean,"  J.  Acoustical  Soc.  Amer.  69  (1)  (1981). 

11.  R.B.  Adler,  L.J.  Chu,  and  R.M.  Fano,  Electromagnetic  Energy  Transmission  and  Radia¬ 
tion,  M.I.T.  Press,  Cambridge,  Massachusetts  (1973). 


12.  L.L.  Beranek,  Acoustics,  McGraw-Hill  Book  Co.,  New  York  (1954). 


-  24  - 


13.  K .  Aki  and  P.G.  Richards,  Quantitative  Seismology  Theory  and  Methods,  W.H-  Freeman 
and  Co.,  San  Fransisco  (1980). 

14.  Sommerfeld,  Partial  Differential  Equations  in  Physics,  Academic  Press,  New  York  (1964). 


-  75  - 


CHAPTER  II: 

THE  HANKEL  TRANSFORM 


H.l)  Overview 

The  relation  of  the  Hankel  transform  to  the  two  dimensional  Fourier  transform  of  a  circu¬ 
larly  symmetric  function  makes  it  as  important  a  tool  for  problems  cast  in  cylindrical  coordinate 
systems  as  the  Fourier  transform  for  proplems  in  cartesian  systems.  Applications  can  be  found  in 
such  diverse  fields  as  astronomy,  electrodynamics,  electrostatics,  oceanography,  physics,  and 
seismology.  Because  it  relates  the  pressure  field  associated  with  a  point  source  in  a  horizontally 
stratified  medium  to  the  plane  wave  reflection  coefficient,  it  forms  the  foundation  of  this  thesis. 
In  this  chapter  we  explore  the  properties  of  the  Hankel  transform. 

We  begin  by  presenting  the  most  common  definitions  of  the  Hankel  transform  in  Section 
n.2.  We  show  haw  the  Hankel  transform  arises  from  the  two  dimensional  Fourier  transform  of  a 
circularly  symmetric  function  in  Section  n.3.  To  relate  the  Hankel  transform  to  the  more  fami¬ 
liar  one  dimensional  Fourier  transform,  in  Section  n.4  we  present  its  asymptotic  form.  In  Section 
II.S  we  complete  our  presentation  of  available  properties  with  a  summary  of  important  results 
available  in  the  literature. 

The  remainder  of  this  chapter  is  devoted  to  results  previously  unavailable.  We  derive  these 
results  to  provide  the  foundation  for  our  work  later  in  this  thesis.  Section  II. 6  examines  window¬ 
ing  and  the  Hankel  transform.  We  will  later  use  the  results  derived  in  this  section  to  determine 
the  range  over  which  pressure  field  data  must  be  known  in  order  to  successfully  estimate  the 
plane  wave  reflection  coefficient.  We  will  also  later  use  the  approximate  results  presented  in  this 
section  to  determine  the  effect  of  varying  source-height  during  data  acquisition  on  the  estimate  of 
the  plane  wave  reflection  coefficient.  Section  n.7  studies  the  effect  of  sampling  on  the  Hankel 
transform.  Sampling  issues  arise  both  when  data  to  be  transformed  is  available  only  on  a  discrete 
set  of  points  and  when  the  Hankel  transform  is  computed  numerically.  The  results  from  this  sec¬ 
tion  will  be  used  extensively  in  Chapter  TV  ■ 
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The  addition  of  white  Gaussian  noise  is  often  a  reasonable  model  for  the  accumulated 
effect  of  many  sources  of  corruption  acting  on  a  measured  signal.  Section  D.8  discusses  the 
degradation  introduced  into  the  Hanlcel  transform  of  a  signal  by  the  addition  of  white  Gaussian 
noise,  ft  also  shows  that  sampling  such  a  function  on  a  square  root  grid  can  improve  the  noise 
behavior  of  the  associated  Hankel  transform . 

We  begin  now  by  presenting  common  definitions  of  the  Hankel  transform. 


0.2)  Definition  of  the  Hankel  Transform 

In  the  literature  a  number  of  different  integral  transforms  are  referred  to  as  the  Hankel 
transform.1  Three  of  these  are  presented  below: 


1  )HT 
2)HT  j 


ryj0(pr)rdr  3  F2(p)  Watson  [1966] 

r)/o(2irpr)rdr  *  Fj( p)  Bracewell  [1965] 


jf(r)j-2w//( 

3)//T3  (r)J  -  }f(r)J0(pr)V^dr  *  F3(p)  Bateman  [1953] 

Definitions  (1)  and  (2)  are  only  superficially  different  since  F 3(p)  =  2rrF1(2-irp).  Definition  (3) 
is  substantially  different  with 


Fj(p)  m 

-TTT-HTx 


(4) 


As  we  will  see,  under  definition  (3)  the  Hankel  transform  has  properties  very  similar  to  the 
Fourier  transform.  We  will  use  definition  (1),  never-the-less,  because  of  its  relationship  to  the 
two  dimensional  Fourier  transform. 


1)  Sometimes  these  trenrfonni  are  also  referred  to  ai  the  zero-order  Hankel  transform.  We  will  not  make 
Mat  distinction  in  this  thesis. 
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U.3)  The  Hankel  Traiuform  as  a  Two  Dimensional  Fourier  Transform 
If  we  use  the  definition  of  Watson 

(')  f  *  ff(r)Jd(Pr)rdr  3  f(p)  (J) 

)  o 

then  the  Hank  el  transform  is  amply  related  to  the  2  dimensional  Fourier  transform  of  a  drcu- 

I 

larly  symmetric  function.  {1,2  ]  To  Aow  this  we  write  the  2  dimensional  Fourier  transform  in 
cartesian  coordinates: 

Pc (kx  /  / fc  (X  ,y  )e l(k‘x  +k>y)dxdy  (2) 

Z1T  —30  — X 

If  /c  (x  ,y  )  is  circularly  symmetric  we  can  unambiguously  define 

fir)  "  fc(x,y)  where  r  ■  V*  Vy2  (3) 

Writing  (2)  in  polar  coordinates  we  have: 

MP.4>)  =  (*-*>rdrde  (4) 

^00 

A  change  of  variables  £  m  0 — <t>  shows  that: 

I  *1* 

*r(p,4>)  =  ^///(Oe'-^rdrdS  (5) 

air  o  o 

so  />(p,<b)  is  not  a  function  of  <J>.  We  suppress  4>,  drop  the  subscript,  P ,  and  perform  the  £ 
integration  using 

2v 

*  /o(x)  (6) 

4ir  o 

to  see  that  any  radial  slice  of  the  two  dimension  Fourier  transform  of  the  circularly  symmetric 
function  fc(x,y)  i*  given  by: 

F(p)  «  ff(r)JQ(pr)rdr  (7) 

0 

which  is  the  Hank  el  transform. 

By  considering  the  Hank  el  transform  as  the  two  dimensional  Fourier  transform  of  a  circu¬ 
larly  symmetric  function  we  can  also  relate  the  Hankel  transform  to  the  Abel  transform.  The 


•r~~ 
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Abel  transform  frequently  arises  in  optics,  seismology  and  other  fields.  In  this  formulation  it  will 
appear  as  the  projection  of  a  two  dimensional  circularly  symmetric  function  onto  its  axis. 

We  begin  by  noting  that  the  slice  of  the  two  dimensional  Fourier  transform  in  polar  form, 
F/>(p,0),  equals  the  slice  of  the  two  dimensional  Fourier  transform  in  cartesian  form,  Fc(p,0), 
since  both  functions  represent  the  same  slice  of  the  two  dimensional  Fourier  transform.  When 
fc(x,y)  is  circularly  symmetric  then  its  transform  in  polar  form,  F^(p,<)>),  is  circularly  sym¬ 
metric  and  equal  to  F  (p),  its  Hankel  transform,  as  we  have  shown.  For  this  case  we  can  there¬ 
for  write: 


f(p)  =  Mp,0)  =  Fc(P’°)  =  T"/  ffc(x*y )eipxdxdy 


If  we  perform  the  y  integration  first  we  have 


f(p)  =  2 ffcOuto  (9) 

The  integral  in  y  generates  the  projection  of  fc  ( x  ,y  )  onto  the  x  axis.  We  define  this  projection 
to  be  p  (x  ).  If  we  use  die  circular  symmetry  of  fc  (x  ,y  )  we  can  rewrite  this  projection  as: 

X  X  _  X 

p(x)  =  / fc (x  ,y  )dy  =  f  fcO/'x2+y2,0)dy  =  2  J/c C^x2jry2 ,Q)dy  (10) 

-X  —  X  O 

Or  in  the  cylindrical  coordinate  system 


P(x)  =  2 


mV,1-*2 


Equation  (11)  is  the  Abd  transform  of  /  (r  ).  The  Abel  transform  can  therefor  be  considered  as 
the  projection  of  a  circularly  symmetric  fc  (x  ,y )  onto  the  x  axis.  Since  the  Hanlcd  transform 
was  shown  to  be  the  Fourier  transform  of  the  projection  we  see  that 


F(p)  -tfT^(r)J  =  FT^(r)J  J  (12) 

This  relationship  was  presented  by  Bracewell.  [3  ]  Implementation  of  the  Hankel  transform 
through  Equation  (12)  is  equivalent  to  the  projection-slice  method  proposed  by  Oppenheim, 
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Frisk,  and  Martinez.  [4  J 

Equation  (12)  relates  the  Hankel  transform  and  the  one  dimensional  Fourier  transform 
through  the  Abel  transform.  When  we  consider  only  large  values  of  p,  the  transform  variable, 
an  approximate  relationship  between  the  Hankel  transform  and  the  one  dimensonal  Fourier 
transform  can  be  developed  that  does  not  involve  the  Abel  transform.  This  can  be  done  through 
the  asymptotic  form  of  the  Hankel  transform,  which  we  present  in  the  next  section. 


11.4)  The  Asymptotic  Form 

If  the  Hankel  transform  is  not  dominated  for  all  values  of  p  by  the  behavior  of  the  kernel 
near  the  origin  (as  would  be  the  case  for  8(r)/r  for  example  )  then  the  asymptotic  behavior  of 
the  transform  can  be  studied  by  substituting  in  the  asymptotic  form  for  the  Bessel  function.  If 
we  use  the  asymptotic  form  for  the  Bessel  function  presented  by  Lipschitz  [5, 6  ]  1 


where  1 0 1  2S  1  and  we  keep  only  the  leading  terms  in  x ,  the  Hankel  transform  becomes: 

VTp> (p)  ~  TT'// O’ )c°* (IPr  I- T V^rdr  (2) 

vrt  o  ** 

If  we  expand  the  cosine  term  Equation  (2)  becomes: 


Vjpj^  (P)  ~ 


_  1 


V7 


//  (r  )cos  (pr  Y^rdr  +  sgn  (p)  //  (r  )sin(pr)  '^r  dr 
o  0 


(3) 


The  integrals  in  Equation  (3)  are  the  Fourier  cosine  transform  and  the  Fourier  sine  transform  [8 
].  hi  some  cases  this  form  allows  us  to  extend  results  available  for  these  Fourier  transforms  to 
die  Hankel  transform.  When  the  sgn  (p)  term  can  be  ignored,  for  example,  Equation  (3)  sug¬ 
gests  that  asymptotically  p  |F  (p)  behaves  much  like  the  Fourier  transform  of  ^  \r\f  (r). 
The  sgn  (p)  term  can  not  be  ignored  without  further  approximation  when  the  values  of  the  sine 


1)  A  n>or«  recast  reference  in  t be  form  of  an  asymptotic  aeries  with  the  tame  leading  term  ii  [7  ] 
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and  cosine  transforms  for  negative  p  effect  the  positive  part  of  the  spectrum.  Such  is  the  case 
when  these  transforms  are  degraded  by  sampling  or  integration  to  a  finite  limit,  for  example.  {9  ] 
Had  we  used  the  definition  of  Bateman  for  the  Hankel  transform,  Equation  (3)  would  have 
appeared  even  more  like  a  Fourier  transform: 

f3(p)  ~  77“ 

v  IT 

Bateman’s  definition  (Equation  II. 2. 3)  is  more  directly  related  to  the  Fourier  transform  than  the 
definition  of  Watson  (Equation  n.2.2).  Despite  this,  we  use  the  definition  of  Watson  became 
we  with  to  preserve  the  relationship  between  the  Hankel  transform  and  the  2-dimensional 
Fourier  transform  presented  in  Section  II. 3. 


J/(r)cos(pr)dr  +  sgn  (p)//(r)sin(pr)dr 
o  0 


(4) 


I 

! 

t 


I 


i 

i 

\ 


1 
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n.S)  General  Properties  of  the  Hankel  Transform 

A  number  of  properties  for  the  Hankel  transform  are  readily  available  in  the  literature. 
[10, 1, 6, 7]  We  present  some  of  the  more  important  of  these  here  for  completeness.1 


PROPERTY 

Hr)  F(p) 

30 

=  jf(rVo(f>r)rdr 

0 

self  -inverse 

F(fi) 

Hr) 

linearity 

o  fi(r)+f2(r)  a 

F  i(p)+F  2(p) 

scaling 

f(ar ) 

-V<*> 

a  a 

derivative 

V2f(r) 

-p  2f  (p) 

power 

X  X 

$f(')gm(')rdr  =  J>(p)G*(p)pdp 

0  0 

moment 

,W-.V 

X 

J r*f(r)dr 

0 

In  the  remainder  of  this  chapter  we  develop  properties  of  the  Hankel  transform  not  avail¬ 
able  in  the  literature  but  of  considerable  importance  to  the  later  developments  in  this  thesis.  We 
begin  by  determining  the  effect  of  on  the  Hankel  transform  of  a  function  when  it  is  multiplied  by 
a  range  limited  window. 

n.6)  Windowing  and  the  Hankel  Transform 
a)  An  exact  windowing  expression 

The  definition  of  the  Hankel  transform  has  infinity  as  the  upper  limit  of  integration.  In 
practice  it  is  often  impossible  to  carry  out  the  integration  to  infinity.  This  may  be  because  the 
function  to  be  transformed  is  only  known  out  to  a  finite  range  or  because  the  integration  must  be 


^We  will  consider  two  functions  to  be  equal  if  the  result  of  oonvoWing  their  difference  with  a  band-ttmitad 
function  is  always  aero.  This  is  equality  in  the  tense  of  generalised  functions. 
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performed  numerically  and  only  a  finite  number  of  calculations  can  be  made.  Following  the 
convention  used  with  the  Fourier  transform  we  will  write  the  upper  limit  of  integration  as  infin¬ 
ity  but  will  make  the  function  to  be  transformed  zero  beyond  some  finite  upper  limit  by  multi¬ 
plying  by  a  window  of  finite  eitent.  [9  ]  In  this  section  we  will  explore  the  degradation  intro¬ 
duced  into  the  Hankel  transform  of  a  function  by  such  windowing.  The  results  of  this  section 
will  also  find  application  to  the  approximate  evaluation  of  integrals  of  the  form 

a 

J/('-Vo(p'-yo«rWr  (1) 

0 

which  arise  in  this  thesis  in  connection  with  source-height  effects. 

An  exact  but  cumbersome  expression  for  the  effect  of  windowing  can  be  derived  from  a 
result  presented  by  Bracewell.  [10  ]  If  we  define: 


J*(P)  *  /p(fy  <jLpr)rdr 
0 

-  (2) 

W(p)  *  jw(r)J^pr)rdr 
0 

Then  Pw  (p),  the  Hankel  transform  of  the  product  of  p  (r )  and  tv  (r )  is  given  by: 


■  x2n  __________ 

*\(p)  m  fp(r)w  (r)Jfor)rdr  =  ffP  (OW  (VV2+£2-2p£cose)Sd0d  §  (3) 

0  0  0 

We  can  relate  the  Hankel  transform  of  the  windowed  function,  Pw  (p)  to  the  Hankel  transform 
of  the  unwindowed  function,  P  (p),  by  carrying  out  the  theta  integration  in  Equation  (3)  to 
obtain: 


*  2v 

J*„(p)  =  (p,4)<*$  with  h( p,0  *  sjV(Vp2+^-2p$cose)</e  (4) 

0  0 

If  H(p,5)  had  the  form  H  (p~5),  then  Equation  (4)  would  be  a  convolution,  reminiscent 
of  the  windowing  result  for  the  Fourier  transform.  [9  ]  By  placing  some  restrictions  of  tv  (r)  we 
can  derive  an  approximate  expression  for  the  effect  of  windowing  that  has  the  form  of  a  convo¬ 
lution. 
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b)  Approximation  as  a  convolution 

A  simpler  approximate  expression  describing  the  effect  of  windowing  can  be  derived  by 
using  the  asymptotic  expression  for  the  Hankel  transform: 


vq^Xp)  ~  —  \ff(r)cos(pr)^rdr  +  sgn  (p) //  (r )sin (pr y^rdr 

w  [o  0  I 

■  [c,(p)  +  Jgn(p)G<j(p)  ] 

G,  (p)  and  Gq(p)  are  a  Fourier  cosine  transform  and  Fourier  sine  transform  respectively.  [8  ] 
The  cosine  transform  and  the  sin  transform  each  have  the  property  that  the  transform  of  a  pro¬ 
duct  is  the  convolution  of  the  transforms.  Using  this,  the  effect  of  windowing  in  the  asymptotic 
formulation  of  Equation  (1): 


V'ipj>„ (p)  ~  —  J/(r )w ( r  )cos (pr  ) v7cfr  +  sgn (p )ff  (r  )w (r )sin  (pr )^7dr  (2) 

tr  o  o 


can  be  written  as: 


VRFw(p)  =  G,(p)*WF(p)  +  sgn (p)  [g„(p)*vMp)J 


where  we  have  defined: 


WF(p)  -  fw(r)e^dr  (4) 

In  general  Equation  (3)  can  not  be  rewritten  as  the  convolution  of  F  (p)  with  a  window 
term  because  of  the  sgn  (p)  term.  However,  if  the  Fourier  transform  of  the  window,  Wp( p),  is 
effectively  confined  to  a  narrow  band  around  p  =  0,  then  for  p  larger  than  this  band  (Bw): 

p  >BW  Ge  (p)* WV  (p)  +sgn  (p)  [g0(P)*HV(p)]  «  Ge(p)*WF(p)+  [sgn  (p)G„ (p)  ]*W>(p05) 
Combining  Equations  (2),  (3),  and  (S)  we  have: 


vTpI>„(p)  *  [vTp]>(p)]*wV(p) 


which  is  our  asymptotic  result. 

If  VFf(p)  **  not  negligible  beyond  some  band,  Bw ,  then  the  effect  of  windowing  on  the 


i 


HanJcel  transform  can  not  be  put  in  this  simple  form.  For  such  cases  the  exact  result  of  Equation 
(II. 6a. 4)  must  be  used, 
c)  Resolution  and  leakage 

Given  Equation  (II. 6b. 6)  we  can  address  the  practical  issues  associated  with  windowing. 
As  is  frequently  done  for  the  Fourier  transform,  we  divide  the  issues  associated  with  windowing 
into  two  general  classes.  The  first  we  call  resolution  and  refers  to  the  local  smearing  affected  by 
the  main  lobe  of  the  window.  The  second  we  call  leakage  and  refers  to  the  contribution  of  the 
side  lobes.  [9, 11  ] 

We  begin  by  expanding  Equation  (I.6b.6)  to  write: 

3D  3D 

p  >0  vftr^p)  =  JV$F(*)Wf  (p-*)d*  With  HV(p)  *  fw(r)e^dr  (1) 

0  -* 

When  p  is  sufficiently  large  (  p  greater  than  some  pq)  then  can  be  considered  constant  over 
the  main  lobe  of  WF  (p-£).  For  these  p.  Equation  (1)  can  be  written  approximately  as 

31 

p>po  (2) 

0 

to  that 

3D 

J’wW®  J>(€)W>(P-«V«  (3) 

0 

Under  this  condition  the  issues  of  resolution  for  die  Hankel  transform  are  the  same  as  those  for 
the  Fourier  transform.  If  we  desire  to  resolve  events  in  the  Hankel  transform  on  the  order  of  5 
then  the  lobe  of  our  window  must  be  less  than  5.  Discussions  about  a  variety  of  windows  are 

available  in  the  literature.  [9, 11  ]  For  the  Hanning  dass  of  windows,  the  main  lobe  width  is 

B 

where  B  is  the  length  of  the  window.  Our  requirement  for  resolution  of  events  on  the  order  5 
becomes: 

* >  ?  w 

Leakage  is  the  phenomena  we  aMOciate  with  the  side  lobes.  For  the  purpose  of  this  analysis 
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WC  ^  ^  **“  £°  ^  ^  *  «  »e  approximated  by  „  topulse  M 

*“  "  ‘8n°re  ^  ^  *“  -  to  resolution.  We  approbate  W,M 

as  a  weighted  superposition  of  impulses: 


wr(  p)  =*  2<j/S(p-r,.) 


The  a,  indicate  the  rate  at  which  the  side  lobes  approach 


becomes: 


zero.  The  convolution  of  Equation  (2) 


30 

^(P)  -  /VtF^Sa^p-S-r,)^ 

-*  l 

^V^FCp-T,)  (6) 

Whenwe^are  concerned  about  the  lealtage  due  to  a  singularity,  Wc  must  ^  ^ 

^  ^  ^  ^  affl°UOt  °f  ,Mka*C  °f  80  -ent  at  f,  of  strength  1  wouid  have  at 

P.  He*  to  Henkel  tomtom  toto ,  *.  p„*,  ^  v_  ^ 

,toW‘  <»  «•"—  leakage.  to  loto  „  ^  „„ 

“  *  *“  ^  to  *«  -1—  to  equivalent  .natonmnee  in  to  p^ 

transform.  For  this  reason  we  have  concentrated  on  the  Hanning  window  rather  than  the  Ham- 
tning  in  many  of  our  examples. 

By  weighting  the  side  lobe  heights  by  a  factor  of  VT  J 

oy  i actor  of  vp>  optimal  wlndows  u)uId  be  desigllcd 

lo,  Henkel  wte.  i,  .  p^ 

d)  Examples 

la  to  .too.  .e  Ptoto  „„  examples  stowing  „„  to  Hmka  amSom  T„  ^ 
we  apply  a  rectangular  window: 

w(r)  _  f1  0  <  '  <  4000 

~  (0  4000  <  r  (1) 

.Pick  to  .  leng*  ,*»„  „  „  t„se  over  sstoh  g.n,  i,  ^ 

m  Sectton  (1.6).  The  first  function  we  transform  is  ~r~  for  which  the  true  Hankel 
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transform  is  given  by  ■  ■  el  **  Figure  IL6d.l  presents  the  log-magnitude  of 

^  p 

eitVr2+(2)J 

— T  '  0  <  r  <  4000.  As  can  be  seen  in  the  figure,  this  function  decays  almost  four 

Vr2+(2)2 

orders  of  magnitude  over  the  window  length.  Figures  II.6d.2a  and  U.6d.2b  present  the  magni¬ 
tude  and  phase  of  its  computed  transform.  Essentially  no  degradation  due  to  aliasing  is  apparent 

ft»VrJ+(  t33)J 

in  the  computed  transform.  Figure  II.  6. 3  presents  the  magnitude  of  -  -  ■  -  ».  Over  the 
window  length  this  function  has  decayed  roughly  two  orders  of  magnitude.  Figures  n.6d.4a  and 


II.6d.4b  present  the  magnitude  and  phase  of  its  Hankel  transform.  The  correct  transform  is 
given  by:  ~  “ 2e>  *2_P^U3^.  The  magnitude  of  the  correct  transform  should  look  like 


for  0:Sp:£*.  Instead  the  magnitude  of  the  transform  shown  in  Figure  n.6d.4a 


shows  considerably  more  degradation  than  that  of  Figure  II.6d.2a.  This  is  due  to  the  fact  that 


this  is  the  transform  of  a  function  which  has  proportionally  more  energy  outside  the  window. 


One  is  tempted  to  assume  the  ripples  apparent  in  Figure  n.6d.4a  arc  due  to  leakage  of  the 
— y  -  —  ■  singularity.  This  is  not  the  source  of  degradation,  however  as  may  be  seen  by  noting 

V*2-p2 

that  this  same  singularity  is  present  in  the  first  transform  of  Figure  n.6d.2a  for  which  no  such 
rippling  is  apparent.  The  rippling  is  due  to  the  smearing  of  the  transform  in  Figure  n.6d.3a  over 
its  rapidly  oscillating  phase  term  e,V^*,_p2^133^  which  is  not  apparent  in  the  magnitude  plot. 
When  the  true  phase  varies  rapidly  over  the  width  of  the  main  lobe  of  the  window  the  effect  of 
smoothing  can  actually  be  to  introduce  rippling  into  the  computed  magnitude. 
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II. 7)  Sampling  anti  Aliasing 

It  is  often  necessary  to  approximate  the  integral  in  the  Hankel  transform  by  a  sum.  This 
approximation  may  be  necessary  because  the  function  to  be  transformed  is  known  only  on  a 
discrete  set  of  points  or  because  the  integral  must  be  evaluated  numerically.  The  resulting  sum 
will  be  a  degraded  version  of  the  true  Hankel  transform.  We  will  adopt  the  terminology  of 
Fourier  transforms  and  refer  to  the  replacement  of  the  integral  by  a  sum  as  sampling  and  the 
resulting  degradation  as  aliasing.  In  this  section  we  examine  the  form  that  aliasing  takes  for  the 
Hankel  transform. 


The  discrete  sum  approximation  that  we  will  concentrate  on  is  the  Fourier-Bessel  series. 
We  will  derive  an  expression  that  relates  the  output  of  the  Fourier-Bessel  series  to  the  true 
Hankel  transform.  Because  the  Fourier-Bessel  series  uses  samples  on  a  set  of  points  that  is 
approximately  evenly  spaced,  the  results  we  derive  will  be  approximately  correct  for  any  evenly 
spaced  sampling  scheme. 

We  begin  with  the  formulation  of  the  Fourier-Bessel  series  [7, 6  )  which  states:1 


. ^  , 

0<P<1  F(p)  =  2XJ - ,-2-  : - JoiKp) 

«  =  io  J  i(X„) 

Where  X.  n  =1,2,3,  •  •  ■  are  the  ordered  zeros  of  Ja(*)- 


If  F  (p)  =  0  for  p  >  1  then  the  integral  in  the  expression  above  is  just  the  Hankel 
transform  of  F(p)  evaluated  at  X„ ,  /(\„)  so  that  the  Hankel  transform,  F  (p),  can  be  expressed 
exactly  as  a  sum: 


TO  f  (\  ) 

0<p<l  F(p)  =22  ~J7~J o(K* P)  when  F( p)  =  0  for  p  >1  (2) 

n-\J  l(*«) 

When  F  (p)  is  not  truly  bandlimited  to  p  <  1  and/or  the  sum  is  not  carried  out  to  infinity.  Equa¬ 
tion  (2)  is  only  an  approximation  to  the  Hankel  transform.  The  study  of  the  effect  of  finite  N  on 


1)  Wa  will  call  two  functions  equal  it  the  Fourier  tnnsfonn  of  their  difference  hat  no  cnerjy  at  any  finite 
frequency.  For  this  reason  we  need  not  sinfle  out  the  values  of  F  (p)  in  Equation  (1)  at  points  of  discon¬ 
tinuity. 
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the  approximation  is  the  study  of  windowing,  covered  in  the  previous  section.  Here  we  consider 
only  the  degradation  that  occurs  because  the  infinite  series  is  used  in  place  of  the  integral. 

Finally,  wt  note  that  it  is  because  the  zeros  of  Jq(x),  A„ ,  rapidly  approach  it  it +-7  that  the 

4 

sampling  above  is  approximately  evenly  spaced. 

To  determine  the  effect  of  approximating  the  Hankel  transform: 


f(p)  =  //(rVo(prWr 
0 


by  the  Fourier-Bessel  series: 


0  <  p  <  1  F(p)  =  £  72^— •/(X.y0(A,p)  (4) 

n  =  lJ  1  (A,) 

we  express  F(p)  in  terms  of  the  correct  transform,  F  (p),  by  inverting  (3)  to  write /(r)  in  terms 
of  F  (p).  We  substitute  this  into  Equation  (4)  to  yield: 


0  <  P  <  1  F( P)  =  2  73~-  J>(5yo(M)&H  'o(*.p)  (5) 

■  [0 

Interchanging  the  order  of  integration  and  summation  we  have 


f  (p)  =  J>«)Mp.ow« 

0 

Where,  following  the  notation  of  Watson  (page  S82)  [7  ]  we  define: 

_  ,  _ £  f  JdtKiVdKp) 

tn  (p4)  -  2  2  - 757— — 

JiKK) 

The  study  of  aliasing  for  the  Fourier-Bessel  series  is  the  study  of  Tx( p,£).  We 
an  expression  for  r*(p,|)  by  using  an  asymptotic  result  presented  by  Schlxfli:1,2  [12  ] 

*r  /_  r\  1  f»hi4»(p-f)  sin  A*(2-p-{)  ]  L  .  _ 1 „ 

Ts (p,5)  “  ctt  — - - W/V_„_n  wh“*  *"  (w 

2VP«  i,i(p-{)  rin  *P.-P-0  4 


1)  This  diflm  from  Watson's  presentation  of  Schlafli'!  result 

2)  SdilafH  does  not  restrict  the  region  of  validity  of  his  result  Watson,  however,  states  that  Schlsfh'i  result 
proceeds  bom  a  formula  which  is  strictly  valid  only  for  0<p+f<2  and  p#{.  We  will  liter  plot  T^(p,( ) 
to  show  that  the  results  of  this  analysis  appear  valid  inside  die  region  0<p+£s2  end  epproahnstety  valid 
outside  that  region. 


I 


ll 


1 

> 


As  N  Tn( p,4)  approaches  a  weighted  sequence  of  impulses.  We  determine  that  sequence 
here. 

We  begin  our  analysis  of  Equation  (8)  by  first  considering  the  expression: 


sin 

sin  Anx 

\r  it* 

N  TTX  +  — * 

4 

««  mx 
sin- 

7rjr 
sin — 

2 

1 

Which  equals 

tin  Nnx  .  itx  , 

- cos  — —  + 

itx  4 

sin  — — 

2 

cos  N  itx 
——  ~  ■  sin 

sin  — — 

2 

In  Appendix  I  we  show  that  as  W  -  *  the  first  term  in  (10)  approaches  the  limit: 


2(-l)*8(f  -2k) 

k  z 


(9) 


(10) 


(11) 


La  Appendix  I  we  also  show  that  the  second  term  in  Equation  (10)  approaches  0.  The  limit 
of  Equation  (9)  is  therefore  given  by: 

Mm  =  2(-D*8(f -2 *)  (12) 

w  -*  nx  *  2 

sin——  * 

2 

Using  Equation  (12)  the  first  term  in  Equation  (8)  can  now  be  seen  to  approach  the  limit: 

_sinAN(p-4)  _  22(-1)*8(p-«-«)  (13) 

sin  Mg-t)  * 

2 

The  second  term  in  Equation  (8)  can  be  put  in  the  form  of  Equation  (9)  by  defining 
y  -  2-p-fc 

sin  An( 2-p-g)  _  sin  ANy 


sin  2LC?-J>-4)  ^  SL 

2  2 


(14) 


Combining  Equations  (13)  and  (14)  we  have: 

^  sin  Af,(2-p-£) 


sin  * 

2 


=  22(-l)*5(2-p-i;-4Jt) 


(15) 


T“" 
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We  can  determine  7,(p,£)  by  combining  (13)  and  (15):1 

Km  TN(P,i)  =  [fi(p-€-4Jfc)  -  8(2— p—?— 4*)]  (16) 

If  our  transform  is  not  severely  aliased  so  that  F  (p)  is  negligible  for  p  >  2  then  substituting 
Equation  (16)  into  Equation  (6)  shows  that: 

2 

0<P  <  1  F( p)  =*  }f (4)^-  [s(p-C)  -  5(2-p-0  (17) 

which  equals  for  0<p<l  :2 

^(P)  =F(P)  ~  ^E^(2-p)  (18) 

VP 

We  observe  that  the  aliasing  result  most  directly  relates  y^pF  (p)  to  VP*‘  (p)- 
An  example  of  aliasing  is  presented  in  Figure  II. 7.1  where  we  see  4"''/p  times  the  Hankel 
transform  of  e  generated  with  the  Fourier-Bessel  series.  The  figure  displays  the  aliasing 
terms  generated  by  the  impulses  in  Equation  (16).  In  the  region  0<  p<  2  the  figure  matches 
the  result  indicated  in  Equation  (17)  very  well.  In  the  region  0<  p<  4  the  figure  does  not 
correspond  exactly  to  what  would  be  determined  by  substitution  Equation  (16)  into  Equation  (6) 
indicating  the  limited  validity  of  Schlafli's  result. 

Figure  n.7.2  shows  a  plot  of  2A//pfr ua(p,£)  0<  p<  10  0<  £<  10.  This  picture  sup¬ 
ports  the  accuracy  of  Equation  (16)  for  7x(p,0  for  0<  p+4s  2  and  suggests  that  Equation 
(16)  is  at  least  approximately  correct  over  the  range  of  p  and  ij  shown  in  the  figure. 

We  conclude  this  section  with  a  final  example  of  aliasing  for  the  Hankel  transform  that  will 
play  an  important  role  in  the  generation  of  synthetic  data.  Figure  II. 7. 3  shows  the  function 

fir)  =  2~"~2  M'o)  >  0  (**) 

r  —r0 

corresponding  to  two  poles,  one  at  r  =  r0  and  the  other  at  r  =  —r0.  This  function  has  the 

known  Hankel  transform:  {8  ] 

1)  Om  thi  ration  of  validity  for  ScUafU'i  result. 

2)  We  have  included  the  point  p  +  (  =  2,  which  if  not  itrictly  within  the  interval  ipedfied  by  Watson. 
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*(p)  =  -ytf<$VoP) 


(20) 


Asymptotically  H^(r „p)' 


V2 


-  i,i 

*  _*>oP 


virr 


'0P 


so  that  its  magnitude  should  appear 


V2  1 

,  Consequently  the  magnitude  of  F  (p)  should  appear  smooth  and  decay  as  — In 

/  WoP  V  D 


p 

Rapid 


Figure  II. 7. 4  we  see  the  numerically  computed  transform  using  the  samples  / 


2048 


oscillations  are  apparent  which  are  not  present  in  the  magnitude  of  the  correct  transform.  The 
source  of  these  oscillations  is  aliasing,  which  can  be  seen  by  using  Equation  (18)  to  approximate 
the  numerically  computed  transform: 


* (p)  -  (2,) 


—iv 


-rrr0p 


Vp  Vrrr0(4096-p)<f 


which  equals 


with  A  s  e(409<5)i,«  jjjjj  cjj,  |,e  rewritten  as: 


^(p)  ~  y  4  +  «n(rop)] 


(22) 


(23) 


'  p  Virro 

where  the  beating  caused  by  the  aliasing  is  apparent  in  the  sin(rgp)  term.  The  aliased  output 

displays  the  form  of  the  —7-  decay  term  times  an  extremely  degraded  estimate  of  >^pF(p). 
Vp  ' 

When  the  transform  decays  only  at  a  rate  of  of  — this  example  shows  that  severe  degradation 

VP 


due  to  aliasing  can  be  expected. 


fl.8)  The  Effect  of  Additive  White  Gaussian  Noise  on  the  Hankel  Transform 
a)  Statement  of  the  Effect  of  Additive  White  Gaussian  Noise  on  the  Hankel  Tranfform 

In  practice  it  is  seldom  possible  to  know  exactly  the  function  whose  transform  we  desire. 
Frequently  it  is  povible  to  model  the  uncertainty  about  the  input  function  by  assuming  that  the 
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errors  associated  with  each  sample  are  random  and  unco rrela ted  from  point  to  point.  Since  the 
combined  effect  of  many  random  factors  can  often  be  modeled  with  a  Gaussian  distribution  (by 
invoking  the  central  limit  theorem  [13  ])  the  assumption  is  often  added  that  the  distribution  of 
error  around  each  point  is  Gaussian.  This  model  of  the  uncertainty  corresponds  to  additive  white 
(the  uncorrelated  assumption)  Gaussian  noise.  We  assume  the  mean  of  the  noise  process  is  zero 
so  that  the  expectation  of  the  noisy  input  signal  is  the  true  input.  If  we  further  assume  that  the 
variance  of  the  noise  process  is  not  a  function  of  the  input  sample  number  then  this  Gaussian 
noise  process  is  stationary. 

In  this  section  we  explore  the  effect  of  such  uncertainty  on  the  Hankel  transform  of  the 
input  function.  Since  the  Hankel  transform  is  a  linear  operator  and  the  noise  process  has  zero 
mean,  the  effect  of  the  noise  will  be  to  introduce  a  variance  in  the  output  of  the  Hankel 
transform  proportional  to  the  noise  power  but  the  expected  output  will  not  be  corrupted.  [14  ]  In 
this  section  we  first  show  that  unlike  the  Fourier  transform  the  variance  of  the  Hankel  transform 
of  stationary  white  Gaussian  noise  is  not  stationary,  but  instead  concentrates  power  near  the  ori¬ 
gin.  This  result  is  important  because  frequently  the  Hankel  transform  is  used  in  place  of  the  two 
dimensional  Fourier  transform  in  problems  with  an  underlying  circular  symmetry.  Because  of 
this  property,  a  slice  of  the  the  two  dimensional  Fourier  transform  of  noisy  measurements  made 
over  a  two  dimensional  grid  of  a  circularly  symmetric  field  will  differ  from  the  Hankel  transform 
of  a  slice  of  that  field. 

In  Section  (b)  we  will  show  that  if  /  ('^T )  is  a  stationary  white  Gaussian  noise  process  then 
jrpft  will  also  be  stationary  white  Gaussian  noise.  This  result  implies  that  if  the  input  function 
is  sampled  on  a  v7  grid  and  each  sample  is  independently  corrupted  by  (zero  mean)  Gaussian 
noise  that  does  not  depend  on  the  sample  number,  then  samples  of  the  Hankel  transform  on  a 
Vp  grid  will  be  independently  corrupted  by  Gaussian  noise  and  the  amount  of  corruption  will 
not  depend  on  the  value  of  p.  On  these  grids  each  sample  represents  die  same  area  of  the  under¬ 
lying  two  dimensional  circularly  symmetric  function  and  the  noise  properties  of  the  Hankel 
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transform  are  equivalent  to  the  noise  properties  of  the  underlying  two  dimensional  Fourier 
transform. 


To  show  that  the  Hank  el  transform  concentrates  noise  power  near  the  origin  we  first  write: 

^a(p)  ■  /[/(0+"(r)]-'o(PrW»-  (1) 

Where  we  have  introduced  the  limit  of  integration,  B ,  to  insure  convergence.  n(r)  is  stationary 
white  Gaussian  noise  with  variance  Nq.  The  variance  of  Ft  is  given  by 


-  E 
*  i 


[Fi(p))  - 

|F»(p)  -£fo(p))  |J] 
|^‘«(rVo(P',Wr|ij 
[■(“)" (0)  ]jo(pa)J(i(p0)a0^a^0 

•  ■ 

■  /^iVo8(a-3Vo(pa)/o(p3)oPd<«<p 

0 

:  N i  (pa)a2da 


For  p  =  0  Equation  (2)  above  shows  that 


M*[f,(0)]  =  Nofajd  a  =  ~- 


When  p#0 


1  ‘  0  P  0 

In  units  of  normalized  frequency  v  ■  p/B 


(2) 


(3) 


(4) 


VAR  [f,(p)  ]  =  W0poJ (pa)ot:da  =  -jp-f&O  (5) 

which  is  plotted  in  Figure  n.Sa.l.  As  can  be  seen  this  function  decays  rapidly  with  v  so  that  the 
Hank  el  transform  concentrates  noise  power  near  the  origin.  We  can  explain  why  this  noise  pro¬ 
perty  of  the  Hank  el  transform  differs  from  that  for  the  Fourier  transform  by  considering  the 
underlying  two  dimensional  circularly  symmetric  Fourier  transform  represented  by  the  Hank  el 
transform. 
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We  recall  from  Section  (II. 3)  that  the  Hankel  transform  of  a  function,  / (r),  corresponds  to 
a  slice  of  the  two  dimensional  Fourier  transform  of  the  function  /(r,B)  made  by  sweeping  /  (r) 
around  the  origin  in  two  dimensions  (so  that  f(r  ,ft)  =  /(r)  for  all  B).  When  we  generate  the 
Hankel  transform  of  the  noisy  input,  f(r )  +  n(r),  we  obtain  a  slice  of  the  two  dimensional 
Fourier  transform  of  /(r)+n(r)  swept  around  the  origin.  The  result  is  very  different  from 
sweeping  /(r)  around  the  origin  and  then  adding  SWGN  (stationary  white  Gaussian  noise)  in 
two  dimensions.  In  the  first  case  the  noise  field  is  circularly  symmetric,  in  the  second  case  it  is 
not.  It  is  the  symmetry  in  the  underlying  noise  field  implied  by  the  Hankel  transform  that  causes 
the  concentration  of  noise  power  near  the  origin. 

We  will  now  show  that  this  behavior  of  the  Hankel  transform  with  respect  to  noise  can  be 
averted  by  changing  to  a  Vr  coordinate  system  for  the  input  and  a  v'p  coordinate  system  for  the 
output.  Samples  evenly  spaced  in  these  square  root  coordinate  systems  have  the  property  that  the 
distance  between  any  two  samples  always  represents  the  same  area  of  the  underlying  two  dimen¬ 
sional  (circularly  symmetric)  function.  Each  noisy  sample  of  the  function  and  its  Hankel 
transform  represents  the  same  amount  of  area  in  the  underlying  two  dimensional  spaces.  Conse¬ 
quently  the  noise  properties  are  equivalent  to  those  associated  with  samples  evenly  space  on  a 
cartesian  grid  (associated  with  the  two  dimensional  Fourier  transform). 

b)  Proof  that  if  f(*r)  is  stationary  white  Gaussian  noise  then  F(^ p)  will  also  be  stationary  white 
Gaussian  noise,  where  F(p)  is  the  Hankel  transform  of  f  (r) 

The  proof  that  if  /(VT )  is  stationary  white  Gaussian  noise  (SWGN)  then  F(v'p)  is  also 
SWGN,  consists  of  three  parts  and  a  conclusion.  First  we  will  show  that  for  the  integral 

transform  defined  as  F3(p)  "  // (r)J<}(pr)'fi>rdr  that  /(r)  is  SWGN  if  and  only  if  F2( p)  is 

« 

SWGN.  Second  we  will  use  this  to  show  that  Vr/ (r)  is  SWGN  if  and  only  if  Vp F(p)  is  SWGN. 
Finally  we  show  that  if  Vr/ (r)  is  SWGN  then  /  (Vr  )  must  be  SWGN. 

i)  Proof  that  Fj(  p)  is  stationary  white  Gaussian  noise  if  and  only  if  f{r)  is  stationary  white  Gaus- 


'4 


•  «*»**<» 


I 
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sian  noise,  where  F2(p)  is  the  Hankel  transform  of  f  (r)  as  defined  by  Bateman. 

We  begin  by  recalling  that  the  Hankel  transform  defined  by  Bateman  [8  |  and  presented  in 

Section  (II. 2)  is  given  by  F2(p)  =  o(pr)v/prdr .  Because  this  operator  is  linear,  F2(p)  is 

o 

Gaussian.  For  the  same  reason  if  the  mean  of  /(r)  =  0  then  the  mean  of  F2(p)  is  also  0.  We 
need  only  show,  then,  that 


(1) 


Now 


(2) 


E  [^(pl^Pl)]  =  CS(p,-p2) 

[f3(Pi)Fi(p2)  j  = 

J7  (r V o(Pir I'^Pi 7drjf  (?Vo(pif)V  P'W^I 

0  0  J 

X  X 

=  J>o8(r  ~  ^ V  V o(P2?) v  Pip2r  fdrd  £ 

0  0 

m 

m  ^ov/PtP5po(pt?Vo(p:{)^f 
0 

Pi 

=  No8(pi-p2) 

This  proves  that  /(r)  SWGN  implies  that  F2(p)  is  SWGN.  The  converse  is  proven  by  noting 
that  the  above  integral  transform  is  its  own  inverse  so  that  the  same  argument  applies,  replacing 
/(r)  with  F 2(p). 

II)  Proof  that  '^pF  (p)  is  stationary  white  Gaussian  noise  if  and  only  if^'rfir)  is  stationary:  whir e 
Gaussian  noise,  where  F( p)  is  the  Hankel  transform  of  f(r) 

From  (i)  we  know  that  if  Vr  f  (r)  is  SWGN  then 


f('^r/(r))J0(pr)^rdr  =  ^ff(r)Jffpr)rdr  ■  Vpf(p) 


(3) 


is  SWGN.  Since  the  Hankel  transform,  Jf  (r)j0(pr)rdr  is  its  own  inverse  the  converse  follcr.s  by 

o 

starting  with  v7/  (r ). 


i 


*AStW  *e 
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iii)  Proof  that  '*rrf{ r )  is  stationary  white  Gaussian  noise  if  and  only  if  f  (vr  )  is  stationary  white 
Gaussian  noise 

First  we  show  that  if  £  [/(rO/ (r 2)  j  =  8(r,-r2)  then  £  J/(/[r1])/(/[rJ)  j  =  ■  ^ 

This  assumes  that  Z(r)  and  l~\r )  exist  in  the  region  of  interest. 


Proof: 


/(*['])  = 

X 

=  j7am)8('-o^ 


da 


=  f/(a)8(r—  l~l(a))-. - — 

L  /(rl(a)] 


Using  Equation  (4)  we  can  write  £  | /  (/  [^1])/  (Z  [''2])  ]  as: 

e\f  (/  [r,l)/(/  [r,])  ]- 

U.E  F(  l)/(  j)J  . /[rWlFWl -  - 


f  ft(a,-a2) 


»[ri-/~1(ai)]S[ra-f~l(q?)](<ai4ai 


}  8[r,-r\a)Mr2-l-'(a)] 

-»  /’[/''(a)] 

If  we  define  £  *  /—1(«)  so  that  a  =  l  {$)  and  da  =  i(£)d%  then 


AL 


(4) 


(5) 


(6) 


=  /8(r,— 08  (r2-Q4 
8(rj-r2) 

Hr1) 

We  now  show  that  /  f^r  )  SWGN  iff  £  j/  (rj)/  (r2)|  =  ^  ^ — — 

We  apply  the  result  of  the  first  part  of  this  section  to  our  special  case  by  defining  l(r)  —  r1 
so  that  i(r)=2rdr  With  this  definition  we  see  that  if  some  p(r)  is  SWGN  then  p(r2)  will  have  a 

variance  proportional  to  — — — — .  This  implies  that  if  fCA)  is  SWGN  then  /(r)  has  a  vari- 

r\ 


ance  proportional  to  — — ^ — —  .  Finally  we  note  that  if  £  j/(ri)/(rj)|  =  — ^ — —  then 
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E  JVr!/  (ri)Vri/  (ri)  |  =  S(ri~ri)  therefor  if  /(vT)  is  SWGN  then  V r/(r )  is  SWGN  The 
converse  follows  from  the  same  steps  in  reverse. 

iv)  Conclusion 

In  (i)  and  (ii)  we  showed  that  if  'Srf(r)  is  SWGN  then  \,p/r(p)  will  also  be  SWGN  In 

(iii)  we  showed  that  v'r/ (r)  is  SWGN  if  and  only  if /(v>  )  is  SWGN.  The  argument  of  (iii)  and 

(iv)  also  showed  that  v'JJf  (p)  is  SWGN  if  and  only  if  f(v/p)  is  SWGN.  Thus  we  have  shown 
that  if  / (v7)  is  SWGN  then  F (V^)  must  also  be  SWGN. 

n.9)  Summary 

In  this  chapter  we  have  developed  the  properties  of  the  Hanltel  transform  in  general.  By 
themselves  these  properties  might  be  amply  interesting  but  together  with  the  ability  to  numeri¬ 
cally  perform  the  Hankel  transform,  they  become  tools  for  scientific  research.  In  the  next  chapter 
we  will  discuss  methods  for  numerically  performing  the  Hankel  transform. 
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Chapter  III: 

Computing  the  Itankei  Transform 


III.  1)  Overview 

Because  of  the  importance  of  the  Hankel  transform  there  presently  exist  many  algorithms 
for  its  evaluation.  As  yet,  however,  there  is  no  generally  accepted  algorithm  analogous  to  the 
FFT  for  the  Fourier  transform.  This  chapter  presents  a  survey  of  numerical  Hankel  transform 
techniques.  It  does  not  include  all  the  published  algorithms  but  does  represent  the  major  classes. 
We  begin  in  section  III.2  by  discussing  perhaps  the  oldest  and  best  understood  algorithm,  the 
Fourier-Bessel  series,  that  we  used  to  derive  our  aliasing  result  of  Chapter  11.  Section  HI.  3 
presents  the  backsmear  method  and  an  example  of  an  efficient  class  of  realizations.  Section  III.  4 
discusses  the  asymptotic  transform  method  as  it  has  been  proposed  in  the  literature  and  presents 
new  results  that  can  be  used  to  improve  this  method  or  to  asses  the  error  in  the  standard  realiza¬ 
tion  from  the  output  transform  alone.  Section  III.5  discusses  a  common  combined  transform 
method  and  presents  a  caution  about  its  use.  Section  III. 6  presents  a  convolutional  method  fre¬ 
quently  used  for  electromagnetic  problems  which  require  the  transform  of  smooth  functions  of 
limited  extent.  Section  m.7  discusses  the  projection-slice  method.  In  that  section  we  develop  a 
fast  algorithm  for  generating  the  projections  (shown  to  be  an  Abel  transform),  which  itself  has 
wide  applications.  [1  ]  When  it  is  followed  by  an  FFT  the  result  is  an  efficient  Hankel  transform. 
Both  the  Abel  and  Hankel  transform  algorithms  are  illustrated  with  examples. 


III.2)  The  Fourier-Bessel  Series 


Probably  the  first  proposed  method  for  evaluating  the  Hankel  transform  is  the  Fourier- 
Bessel  series,  which  was  discussed  in  the  aliasing  section.  The  Fourier-Bessel  series  is  summar¬ 
ized  by  the  identity:  [2, 3  ] 


0<P<1 

,-io 


(i) 


Where  n  =1,2,3,  •  •  •  axe  the  ordered  zeros  of  Jq(x). 


This  relation  can  be  used  to  obtain  the  Hank  el  transform  of  a  band  limited  function,  /  (r), 
defined  by 


^(P)  s  ff(r)J0(pr)rdr  (2) 

o 

where  F (p)  =  0  for  p>l  by  noting  that  when  F(p)=0  p>l: 

*  1 

HK)  =  fF(p)JdpK)pdp  -  J>(p)J0(pMp<*P  (3) 

0  0 

Substituting  into  Equation  (1)  we  have: 

*  f  (K  ) 

0<P<1  F  (p)  =  2  2  1^Ljq{Kp)  (4) 

M  =  lJ i  (An) 

which  is  in  the  form  of  a  sum  as  we  desired.  If  F(p)  is  bandlimited  to  p<fi ,  instead  of  p<l, 
Equation  (4)  can  be  modified  by  a  change  of  variables.  The  resulting  general  form  states  that 
when  F  (p)  =  0  for  p>fl  then 


0<p<fl 


*"(P) 


=  JL 

B2, 


B  o(-f-P) 


Hk) 


B 


(5) 


The  properties  of  the  Fourier-Bessel  series  have  been  extensively  studied.  [2, 3  ]  As  a  numerical 
algorithm  for  implementing  the  Hankel  transform  it  is  usually  appropriate  only  when  a  few 
values  of  F( p)  are  required  or  when  computation  cost  is  less  important  than  careful  control  over 
the  errors.  This  is  because  the  Fourier-Bessel  series  requires  than  a  new  sum  be  computed  for 
every  value  of  p.  Further,  it  requires  that  the  function  to  be  transformed  be  available  on  the 


nonuniform  grid, 


hjL 

B  ' 


Another  algorithm  that  is  appropriate  when  only  a  few  values  of  F(p)  are  required,  but 
which  accepts  input  values  on  an  even  grid,  is  the  backsnear  method.  We  present  the  backsmeax 
method  in  the  next  section. 
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ni.3)  The  Backsmear  Method 

The  backsmear  method  exploits  the  efficiency  of  the  FFT  to  provide  an  efficient  algorithm 
when  only  a  few  output  values  of  the  transform  are  desired.  This  algorithm  can  be  derived  by 
replacing  the  Bessel  function  with  its  integral  representation  the  Hankel  transform  can  be  written 
as  the  two  dimensional  integral: 

*  2v 

^(P)  s  ff(r)J0(pr)rdr  =  ~~  f 
0  217  o 

X 

If  we  define  G(x)  =  jf  (r  )re irx  dr  and  note  that  COS0  is  even,  the  transform  becomes: 

0 

V 

v  2  . 

F(p)  =  —  /G(pcos8)d8  =  —  f  G(pcos8)+G(-pcos0)  <f»  (2) 

IT  0  '•*’  0  ^  ■* 

G(x)  can  be  evaluated  efficiently  on  an  even  grid  using  an  FFT.  [4  ]  The  integration  called 
for  in  Equation  (2),  however,  must  be  performed  for  every  desired  value  of  p.  Further,  since  the 
integration  of  Equation  (2)  is  in  8,  some  interpolation  scheme  must  be  used  to  generate 
G(pcos8)  on  the  quadrature  points  required  by  the  numerical  integrator. 

The  backsmear  method  is  appropriate  when  the  transform  is  required  at  only  a  few  values 
of  p.  When  this  is  the  case  and  when  G  (x)  is  slowly  varying  so  that  the  numerical  integration 
which  must  be  carried  out  upon  it  can  be  done  efficiently,  then  the  backsmear  method  is  effi¬ 
cient.  This  is  particularly  true  if  the  interpolation  scheme  used  permits  a  quadrature  formula  to 
be  developed  for  the  numerical  integration. 

One  such  development  is  presented  below  for  the  case  of  linear  interpolation  between  the 
points  of  L(x)  m  G(x)  +  G(-x). 

We  seek 


J/(r),r'>c0,<Wr 


d9 


(1) 


2 

*  0 


St 

F(p)  =  —fL(pcosO)JO 


(3) 


tod  have  available  only  L  (*— *)  from  the  output  of  the  FFT.  We  approximate  L(x)  by  linearly 
NT 
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f 
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! 

i 

i 
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interpolating  L  (~j^)  M  that 


L(x)=  2MX> 

u  j‘Q 


L(m)  +  (X~M)(NT)  L{1W)~L{m)  ~m<x<JNT 

0  otherwise 


We  can  integrate  L  (pcos0)  as  follows 


Jl(pcos0)d0  =  /  2  £/(pcos0)<f@  =  2  J*£/(pcos0)d0 

o  0  j-1  ]-l  o 

The  L,  terms  can  be  integrated  term  by  term: 


1  at 

-  l(M + 

“**  err 

where 

.  _!  (cOS_1Jt  \X  |Sl 
068  x  =  \o  |x  |21 

Evaluating  this  integral  and  substituting  into  Equation  (6)  we  have 


h  (?) 


**  1 

■  i1 


F  (— )  = 
''NT' 


IN* 

•  2  co,_,(i‘) '  c“",(i*L)]  “  4(^ 

+  2  [v*3r?  -  v*2-0+i)3]  [i  (^)  - 


AT  )  ^JVt' 


When  k  is  small  or  /  ( - )  is  desired  for  only  a  few  values  of  k ,  this  is  an  efficient  pro- 

NT 

eedure.  When  k  is  large  and  /(-— ■)  is  required  for  many  values  this  is  a  slow  method,  how- 

NT 

ever.  In  the  next  section  we  present  a  fast  algorithm  for  the  evaluation  of  the  Hank  el  transform 
at  large  values  of  the  transform  variable,  p. 


i  ' 
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III. 4)  The  Asymptotic  transform 

In  the  seismic  community  the  Hankel  transform  is  commonly  evaluated  approximately  by 
replacing  the  Bessel  function  with  its  asymptotic  form  so  that  the  resulting  integral  looks  like  a 
Fourier  transform.  [5,6  ]  The  result  will  usually  asymptotically  approach  the  Hankel  transform, 

though  pathological  functions  can  be  found  for  which  this  is  not  the  case  ( ^  for  example). 

The  main  disadvantages  of  this  method  are  that  it  is  almost  never  suitable  for  small  values  of  p 
and  that  the  error  induced  by  the  substitution  is  determined  in  part  by  the  function  being 
transformed.  The  main  advantage  is  that  the  resulting  integral  looks  very  much  like  a  Fourier 
transform  and  can  be  evaluated  efficiently  with  an  FFT. 

An  expression  for  the  Bessel  function,  suitable  for  an  asymptotic  expansion  is  provided  by 
Lipschitz:1  [7  ] 

for  X  >0  Jq(x)  =  (*---)  +  ~sin(x-~)  -  -^-iroifx --J)  +  ~  |  1®I<1  (1) 

In  the  literature  describing  this  technique  [5  }  only  the  leading  term  is  retained,  providing: 

f(p)  ~  ■^^//('•)cOs(pr--2-)v7</r  when  p  >  0  (2) 

This  is  written  as  a  Fourier  transform  by  noting  that 

COs(pr  “  )  =  -^(cos(pr)  +  sin  (pr))  (3) 

hence 


F  (p)  ~ 

vp*jr 


//(r)V7cos(pr)dr  +  ff  (r  J  sin(pr  )dr 
.0  0 


when  p  >  0 


(4) 


This  expression  can  be  evaluated  with  an  FFT. 

We  can  write  the  error  associated  with  the  transform  technique  as 

V5 


®(p)-//(r)r 


\J o(pr  )  ~  T7-~"  c°r  (pr  ~  ~r) 
irpr  4 


dr 


=  jf/(r)ff(pr)dr 


(5) 


1)  A  more  recent  expansion  with  the  lime  leading  behavior  can  be  found  in  Watson.  [2  ] 


(9) 
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and  note  that 


*  X 

=  -ff(r)Vrsi/t(i>r)dr  = 

*»P  0  V  r  0 

=  //(r)v7coi(pr)dr  =  E£(p) 
°P  o  v  r  0 

.  c 

— j/^^-coi(pr)dr  =  -fc  (p) 

op  o  r 

dp  0  r 


f  |»  p  P  ( 

F{p)  =  Mp)+^(p)  + 


We  can  combine  terms  above  by  defining 


fl(p)  -  ?c( P)  +  -P.(p) 


so  Equation  (11)  becomes 


*■ w  -  vfc-  r(p) + +  (i3) 

The  indefinite  integrals  above  can  be  converted  to  definite  integrals  by  integrating  from  0  to 
infinity  and  adding  on  the  values  necessary  to  force  the  resulting  equation  to  match  Equation  (8) 


at  p  =  0: 


+  C,  with  Ct  »  f-^p-dr  (14) 

0  0  Vr 

//fittiyM«  =  /[/f.(Cy«  +  C,ld?  +  C2  with  C2-/^dr  (15) 


Performing  the  integrations  over  the  constants  and  replacing  the  iterated  integral  we  are  left  with 
the  expression 


The  first  term  can  be  recognized  as  the  usual  asymptotic  expression  for  the  Hankel  transform. 
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The  second  and  third  terms  are  corrections  to  this  expression  which  can  in  principle  be  deter¬ 
mined  from  it.  The  final  term  is  the  remaining  error  term  but  which  is  considerably  smaller  than 
the  term  associated  with  the  usual  asymptotic  expression. 

If  it  were  desired  to  calculate  the  second  and  third  terms  directly  as  a  Fourier  transform, 
the  same  procedure  could  be  used  as  will  be  described  in  Section  (III.7)  where  it  is  applied  to  cal¬ 
culating  the  Abel  transform. 

Equation  (16)  can  also  be  interpreted  as  providing  first  and  second  order  error  estimates 
for  the  classical  asymptotic  transform.  These  estimates  allow  the  error  associated  with  the  classi¬ 
cal  method  to  be  estimated  (to  first  and  second  order)  from  a  knowledge  of  only  the  resulting 
transform.  Such  estimates  permit  one  to  interpret  the  result  of  an  asymptotic  transform  with  a 
questionable  appearance. 

III.5)  A  Combined  Transform  Method 

The  Fourier-Bessel  series  and  the  backsmear  methods  both  permit  the  calculation  of  the 
Hankel  transform  on  any  output  grid.  The  computational  cost  of  each  of  these  methods  increases 
linearly  with  the  number  of  points  computed.  The  asymptotic  method  is  fast  and  usually  gen¬ 
erates  good  estimates  of  the  Hankel  transform  when  p  is  large.  A  combined  scheme  is  possible 
which  uses  a  slow  method  such  as  the  Fourier-Bessel  series  or  the  backsmear  to  compute  the 
Hankel  transform  point  by  point  for  low  values  of  p  and  which  switches  to  the  asymptotic 
transform  for  large  values  of  p.  Such  a  method  has  been  proposed  in  the  literature.  [8  ] 

The  main  issue  with  such  a  scheme  is  the  point  at  which  the  algorithm  should  cease  comput¬ 
ing  the  transform  point  by  point  with  the  slow  method  and  begin  to  accept  the  output  of  the 
asymptotic  transform.  At  present  there  is  no  reliable  method  for  doing  this.  It  has  been  sug¬ 
gested  in  the  literature  that  the  transition  be  made  at  the  first  point  for  which  the  slow  algorithm 
produces  a  transform  value  which  matches  the  value  of  the  asymptotic  transform  within  a  speci¬ 
fied  tolerance.  (8  ]  This  scheme  for  switching  to  the  asymptotic  transform  would  be  entirely  ade- 
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quate  if  the  error  made  by  the  asymptotic  algorithm  were  monotonically  decreasing  in  the 
transform  variable,  p.  This  error  is  not  always  monotonically  decreasing,  however,  as  we  will 
illustrate  by  constructing  a  function  for  which  it  is  not.  The  existence  of  such  a  function  indicates 
that  the  switching  method  proposed  in  the  literature  [8  ]  will  not  always  work.  In  fact,  for  the 
function  we  construct,  the  error  made  by  the  asymptotic  method  will  be  zero  at  a  point  we 
specify.  The  switching  scheme  proposed  would  begin  to  accept  the  asymptotic  transform  before 
this  point  even  though  the  error  made  by  the  asymptotic  transform  beyond  this  point  might  be 
greater  than  the  specified  tolerance. 

As  in  Section  (III.4)  we  write  the  error  associated  with  the  asymptotic  transform  as: 


«(p)  =  7  ff(r)H(pr)dr 
P  o 


where  H  (pr)  *  pr  l./0(pr)  “  ~r—~rCOS  ( pr — — )  1  and  was  plotted  in  Figure  m.3.1.  In 
[  vrrpr  4  J 

order  to  construct  a  function,  f(r),  such  that  e(p)  is  not  monotonically  decreasing  we  first 
choose  some  small  8  and  set 


1  0<  r<  8 

8 

(pnr)dr 
Po  o _ 

X 

7 ~fn(par)dr 
Po  s 


8<  r<  « 


For  this  /  (r  )  the  error  made  by  the  asymptotic  transform  is: 

as 

«(p)  =  ^SH(pr)f(r)dr 

»  (3) 

,  *  ~SB(p(/)dr  , 

=  l/ff(pr)dr~ ^ -  1/ff(pr)dr 

Po  1  "  P( 

■~fH(pgr)dr 

Po  8 

When  p  =  po  the  error,  e(po),  equals  zero.  In  general  for  p  >  pg  it  will  not.  For  this  con¬ 
structed  function  the  error  is  not  monotonically  decreasing  and  the  switching  procedure  described 
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in  the  literature  will  not  work. 

Until  an  adequate  switching  method  has  been  found  or  until  the  class  of  functions  for  which 
the  proposed  switching  algorithm  will  work  has  been  well  defined,  there  can  be  no  guarantee 
that  the  combined  algorithm  will  work  properly  and  this  method  should  be  used  with  caution. 

II1.6)  A  Convolutional  Method 

In  this  section  we  describe  a  method  for  computing  the  Hankel  transform  that  puts  the 
Hankel  transform  in  the  form  of  a  convolution  by  transforming  to  an  exponential  grid.  When 
this  grid  does  not  involve  an  extraordinary  number  of  points  to  adequately  represent  the  func¬ 
tion,  the  Hankel  transform  can  be  efficiently  evaluated  with  an  FFT.  We  will  discuss  the  presen¬ 
tation  of  this  method  by  Siegman  [9  ],  though  other  presentations  are  available  in  the  literature. 
[10] 

Siegman  converts  the  Hankel  transform  into  a  convolution  by  sampling  the  function  on  an 
exponential  grid.  He  begins  with  the  Hankel  transform  integral 

30 

^(P)  =  ff(.r)Jo(.Pr)rdr 
0 

After  the  following  definitions: 

PsPo«'  v  =  /«(-£-) 

Po 

r  m  'a**  ?  =  Mr-) 
r0 

Equation  (1)  becomes: 

F(  PoO  = 

30 

=  Sf(i‘o‘tyo(Wo‘v+i)>‘o*2,i<lS 

—  30 

30 

=  roj*  ['2V('-o'?)J^o((Vo«("+iV5 

30 

=  riS  [.-/(ror-.)]/^^. 


(1) 


(2) 
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which  is  the  convolution  of  r$e  2“f  (rge  ")  with  /o(Port*0  and  can  *1C  implemented  with  an 
FFT. 

The  strength  of  this  technique  is  the  efficiency  with  which  the  Fourier  transform  can  be 
implemented.  Its  weakness  stems  from  the  requirement  that  /(r)  be  sampled  evenly  in  e'*.  In 
order  to  obtain  the  sampling  density  necessary  to  represent  a  function  near  the  origin  it  is  likely 
that  such  a  density  of  points  is  necessary  to  represent  the  function  at  larger  ranges  that  the  com¬ 
putational  savings  are  lost.  Also,  the  presence  of  the  gain  factor  e  ~2u  might  be  a  severe  problem 
for  the  region  0<  r  <  1.  It  is  unlikely  that  this  transform  technique  would  work  efficiently  for 
functions  with  more  than  moderate  range-bandwidth  products. 

In  the  next  section  we  present  another  Hankel  transform  algorithm  that  exploits  the  compu¬ 
tational  efficiency  of  the  FFt  through  a  change  of  coordinate  system.  It  requires  only  that  the 
function  be  represented  on  a  square  root  grid,  however. 

m.7)  The  Projection-Slice  Method 
a)  Overview 

In  Section  II.  2  we  related  the  Hankel  transform  to  the  two  dimensional  Fourier  transform 
of  a  circularly  symmetric  function.  We  showed  that  the  Hankel  transform  can  be  obtained  by 
first  forming  the  projection,  or  Abel  transform: 

p(r)=A  f{r)  m  2//(V^V)<fy  =  (1) 

which  is  then  followed  by  a  Fourier  transform: 

f(p)  ■  =  -z-  fp(r)e lfrdr  (2) 

0  -» 

Oppenhesm,  Frisk,  and  Martinez  [11  ]  suggested  that  the  computational  efficiency  of  the  FFT  be 
exploited  by  implementing  the  Hankel  transform  as  a  numerical  projection  followed  by  an  FFT. 
Previously,  however,  the  projections  were  expensive  to  compute,  requiring  on  the  order  of  n2 
operations  and  function  evaluations  or  interpolations.  As  part  of  this  thesis  a  iomputationally 
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effirient  method  for  computing  the  projections  (the  Abel  transform)  was  developed.  [12  ]  When 
followed  by  an  FFT  the  result  is  a  computationally  efficient  Hankel  transform  requiring  on  the 
order  of  N*logN  operations. 

h)  The  Abel  transform 

We  consider  the  Abel  transform  shown  below: 


As  suggested  by  Bracewell  [1,2],  we  write  this  transform  as  a  convolution  by  defining: 

*(0- 


0  r  20 

1 


r<  0 


'  -r 


and 


p(r)  *p(y7)  ra 0 

/»  -/(V7)  raO 

which  leads  to  the  convolution  formula: 

Pir)  =  fir)  *  h{r) 

The  Abel  transform  of /(r),  pit),  is  determined  by: 

Pir)  ~  Pir*) 

Bracewell  [2]  proposed  evaluating  (6)  by  shifts  and  sums. 

Because  the  Fourier  transform  of  lt(r)  is  the  known  analytic  function: 

H  O')  =  for  all  v 

pir)  can  be  determined  in  principle  by  means  of  the  Fourier  transform: 

1+1  1  I2vr» 


Pir)  =  Jf(v)-±±!- 


(1) 


(2) 


(3) 

(4) 


(5) 


(6) 


(7) 


(8) 


where  F(v)  is  the  FT  of /(r).  Unfortunately  the  angularity  at  v  =  0  makes  this  function  dif- 


u 
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ficult  to  represent  in  the  computer  and  is  responsible  for  the  long  tails  of  p  (r  )  which  cause  alias¬ 
ing  problems  when  the  convolution  is  implemented  with  a  Fourier  transform.  We  avoid  these 

difficulties  by  removing  the  — W  factor  from  the  numerical  part  of  the  transform  (8)  in  such  a 

v  v 

way  that  the  remaining  function  is  as  well  behaved  as  F  (v)  within  the  numerical  portion  of  the 
transform.  The  singular  part  of  the  transform  is  performed  analytically  and  added  in  after  the 
numerical  processing. 

To  this  end  the  transform  (8)  is  written 


p(r)  ■^el2vr“dv  +  ■^el2°r'dv 


=  !±i  r/(v)  -  f *  +  j?(0)jflz^#-»-fa-r).rfv  +(9b) 

2  q  Vy  Q  Vv 

> 

+  /  \~yh.e‘*»  dv  +  F(0)  dv 


Where  b  is  a  parameter  analogous  to  the  real  part  of  the  exponential  in  a  Laplace  transform. 
Our  choice  of  b  is  described  later. 

The  integrands  in  the  first  and  third  integrals  of  F.quation  (9h)  do  not  have  singularities  at 
v=  0.  Because  both  the  numerator  and  denominator  of  these  integrands  approach  0  as  v 
approaches  0,  they  can  be  evaluated  by  1'  Hospital's  rule  to  show  that  as  v  approaches  0  they 
approach  F(0). 


Upon  defining 


r/(0) 

*  1  (F (v)  -F(0)e"*l,|(l-^))Av 


v=0 

otherwise 


and  performing  analytically  the  two  integrals  that  do  not  depend  on  F  (v)  we  have  : 


'<r>-  t1  o) 


1  i^ir 


2irir  +2irir 
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L(v)  was  defined  in  (10)  such  that  Z.(0)=F(0)  and  L  (p)~  F  (v)  for  large  v.  The  parame¬ 
ter,  b  ,  was  chosen  so  that  L  (v)  moves  smoothly  between  its  limits.  If  b  were  set  to  0,  L  (v) 
would  have  a  DC  term  that  would  transform  to  an  impulse.  Theoretically  this  would  be  canceled 
by  the  angularities  in  the  portion  of  the  transform  performed  analytically  (see  equation  above) 
but  computational  errors  would  certainly  cause  problems.  If  b  were  infinite,  L  (v)  would  suffer 
from  the  lA^v  singularity  at  the  origin,  b  is  chosen  to  smooth  out  the  singularity  somewhat 
between  these  limits.  We  have  been  using  values  of  b  such  that  e  ~b“  has  decayed  to  e  -1  after 
roughly  6  samples  of  F  (v). 

We  present  three  examples  of  functions  processed  with  the  Abel  transform  algorithm 
described  above. 


Example  1  (a) 


<  1 
=  1 
>  1 


(12) 


The  first  example  is  the  transform  of  the  pillbox  function 

/(O-fe  I' 

v.0  |r 

1024  samples  of  this  function  were  generated  on  a  vn T  grid  with  T=  1/32  and 
transformed.  The  result  is  shown  in  Figure  m.7b.l  as  dots  superimposed  over  a  solid  line  which 
is  the  transform  computed  analytically  (  2^l-r2  ).  The  output  matches  the  analytic  solution 
well. 


| 

i 

i 


> 


/ 


Example  2  (a) 

Fat  the  second  example  we  transformed  the  function 

^i(0  ,  . 

-J-  *(r)  (13) 

where  w  (r  )  is  a  Hanning  window.  2048  samples  on  a  vn7r  grid  with  T=  1/2  were  input  (Figure 
HI. 7b. 2).  Figure  in. 7b. 3  shows  the  output  (dots)  superimposed  over  the  correct  transform  (solid 
line).  The  correct  Abel  transform  was  computed  by  evaluating  the  Fourier-Bessel  series  [11}  to 
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obtain  a  stow  but  accurate  Hank  el  transform  of  the  windowed  input.  The  Hankel  transform  was 
then  inverse  Fourier  transformed  to  generate  the  Abet  transform.  In  the  absence  of  the  window 
the  result  would  have  been  sin(r)/r. 

The  output  is  coincident  with  the  correct  solution. 

Example  3  (a) 

For  the  third  example  we  again  transformed  2048  points  of 

(20) 

on  the  grid  ^ nT  but  now  T  was  chosen  to  be  4.  This  input  is  shown  in  Figure  m.7b.4. 
Increasing  the  sampling  interval  reduced  the  effect  of  windowing  on  the  input  because  a  greater 
range  of  the  function  was  represented  but  it  also  increased  the  sampling  interval  on  the  output. 
Figure  HI. 7b. 5  shows  the  output  (dots)  superimposed  over  the  correct  transform.  Again,  there  is 
no  discernible  error. 

c)  The  Hankel  Transform 

To  complete  the  Hankel  transform  it  is  necessary  to  Fourier  transform  die  projection 
obtained  from  the  Abel  transformer.  Unfortunately  this  is  available  on  a  v  nT  grid  and  not  the 
even  grid  required  by  the  FFT.  To  generate  p  (r  )  on  an  even  grid  it  is  necessary  to  interpolate. 
If  a  simple  interpolation  scheme  is  used,  like  sample  and  hold  or  linear  interpolation,  the  result 
will  be  generated  rapidly  but  may  suffer  some  degradation.  If  a  more  sophisticated  interpolator 
is  used,  better  results  can  be  expected  but  at  the  expense  of  greater  computation  time.  Because 
the  interpolation  is  from  an  uneven  grid  to  an  even  grid  (and  not  the  reverse)  it  is  difficult  to 
characterize  the  error  theoretically  beyond  the  fact  that  the  finer  the  initial  grid  the  better  the 
results.  We  complete  the  Hankel  transform  of  the  examples  presented  above  using  linear  interpo¬ 
lation  to  generate  the  uniform  grid. 
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Example  I  ( h ) 

Figure  m. 7c.  1  shows  the  result  of  using  a  FFT  on  the  linearly  interpolated  projection  generated 
in  Example  1  (a).  The  dots  are  the  calculated  output  and  are  superimposed  over  the  analytic 
solution  (solid  line).  The  agreement  is  excellent.  The  time  required  to  perform  the  total  Hankel 
transform  (1024  input  points  to  1024  output  points),  including  the  required  linear  interpolation, 
was  less  than  31  seconds  on  a  PDP  11-55  with  a  floating  point  processor. 

Example  2  (b) 

Figure  ni.7c.2  shows  the  result  of  Fourier  transforming  the  linearly  interpolated  output  of 
Example  2  (a).  Again  the  dots  represent  the  output  of  the  Abel-Fourier  scheme  and  the  solid 
line  is  the  Hankel  transform  as  computed  by  the  Fourier-Bessel  series.  [11]  The  agreement  is 
excellent. 


Example  3  (b) 

Figure  IH.7C.3  compares  the  result  of  Fourier  Transforming  the  linearly  interpolated  output 
of  Example  3  (a)  (dots)  with  the  correct  transform  (solid).  Significant  distortion  is  apparent  in 
this  transform.  Since  the  output  of  the  Abel  transform  in  example  3  (a)  essentially  equals  the 
output  in  example  2  (a)  (the  correct  projection)  except  for  the  sampling  interval,  we  can  associ¬ 
ate  this  distortion  with  the  linear  interpolation  performed  before  the  FFT. 

d)  Discussion 

We  have  found,  as  indicated  by  the  examples  above,  that  the  Abel  transformer  works  well. 
When  its  output  can  be  successfully  interpolated  and  is  followed  by  a  FFT  the  result  is  a  fast, 
accurate  Hankel  Transform  as  illustrated  by  examples  1  and  2.  As  the  spacing  between  output 
samples  of  the  Abel  transformer  is  increased,  the  suitability  of  a  simple  interpolation  scheme 
becomes  suspect.  Example  3  was  chosen  to  illustrate  the  effect  of  inappropriate  interpolation  on 
the  resulting  Hankel  transform.  At  this  point  it  would  be  prudent  to  determine  the  validity  of  a 
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Hankel  transform  performed  with  this  algorithm  by  comparing  the  output  for  inputs  of  different 
grid  spadngs. 

e)  Summary 

The  procedure  for  performing  the  Hankel  transform  Hf(r)  =  FH  (v)  is  summarized 
below. 

1)  generate  f  (r)  »  fC^r) 

2)  Fourier  transform  to  obtain  F  (v) 

3)  generate L  (v)  ■  F(v)  -  F(0)e'*l“l(l-'4)/V^ 

4)  perform  the  inverse  Fourier  transform  and  add  in  the  analytic  terms: 

Hr)  -  ffr-tW  -  '«»{vrfer  -  Tkz  -  -  ?rb}} 

5)  interpolate  to  an  even  grid  p(r)  *  p  (r2) 

6)  Fourier  transform  to  obtain  the  Hankel  transform 

*'*(*)  =  FT  p(r) 

Each  of  steps  2)  thru  6)  requires  no  more  than  the  order  of  N  iog(N)  operations.  There¬ 
fore  the  total  transform  can  be  accomplished  on  the  order  of  N  log(N)  operations.  Direct  compu¬ 
tation  of  projections  from  the  2  dimensional  circularly  symmetric  function  would  have  required 
at  least  N  function  evaluations  and  N  sums  for  each  of  N  points  before  the  final  FFT,  which 
leads  to  an  algorithm  requiring  on  the  order  of  N1  operations.  Therefore  for  sufficiently  large  N 
this  method  of  cslmlsting  the  projections  can  provide  a  considerable  advantage  in  speed. 

Steps  1)  and  5)  above  indicate  that  in  two  places  interpolations  may  be  required.  In  many 
cases,  however,  the  function  to  be  transformed  can  be  generated  initially  on  a  grid  evenly  spaced 
faV 7.  Further,  /(VnT)  ,  as  required  by  this  algorithm,  has  the  desirable  feature  that  equal 
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areas  of  the  underlying  two  dimensional  function  /(r, B)  are  represented  between  samples.  If 
stationary  white  gaussian  noise  (SWGN)  corrupts  the  measurement  of  / (y'rtT  )  n=  0,1,  '  '  - 
then  the  Hankel  transform  on  a  grid  will  be  corrupted  by  SWGN  (corruption  of  equal  areas 
of  the  underlying  two  dimensional  function  produces  corruption  of  equal  areas  of  the  underlying 
2-D  FT).  This  is  not  true  if  / ( nT )  is  corrupted  by  SWGN. 

To  implement  a  Hankel  transform  using  this  method  it  is  necessary  to  perform  the  interpo¬ 
lations  of  step  5).  Because  of  the  speed  of  the  Abel  transform  portion  of  this  algorithm  we  have 
found  it  sufficient  in  many  cases  to  simply  generate  an  over  sampled  version  of  p  (r)  and  to  use 
linear  interpolation  to  obtain  p  (/■  )  . 

f)  Conclusion 

For  many  applications  this  method  of  calculating  the  Abel  transform  appears  to  permit  the 
efficient  calculation  of  the  Hankel  transform  for  large  data  files.  This  transform  method  is  par¬ 
ticularly  appropriate  for  the  evaluation  of  the  Sommerfeld  integral,  in  which  the  oscillations  of 
the  kernel  increase  with  the  independent  variable.  As  a  general  transform  method  issues  of 
representation  on  a  grid  must  be  further  explored.  Because  of  the  equal  area  property 
described  earlier  for  / (Vr  )  and  because  the  speed  of  this  algorithm  permits  oversampling  in 
p(r2)  it  is  not  expected  that  these  issues  will  pose  serious  problems.  It  appears  that  the  Vr  grid 
is  of  fundamental  importance  in  the  Hankel  transform 
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Ill.  8)  Summary 

In  this  chapter  we  presented  a  number  of  numerical  techniques  for  evaluating  the  Hankel 
transform.  No  one  technique  is  ideal  for  all  situations.  When  the  value  of  the  transform  is 
desired  at  only  a  few  points,  the  Fourier-Bessel  series  or  the  backsmear  method  is  appropriate.  If 
speed  is  extremely  important,  accuracy  is  not,  and  the  transform  is  not  needed  for  small  argu¬ 
ments,  then  the  asymptotic  method  is  justified.  If  the  input  function  and  its  Abel  transform  can 
be  well  represented  on  a  square  root  grid  (which  is  the  case  for  functions  which  increase  in 


complexity  with  range)  then  the  Hankel- Abe)  (or  projection-slice)  method  is  a  good  choice. 

Having  established  the  properties  of  the  Hankel  transform  and  examined  its  numerical 

implementation  we  are  now  ready  to  consider  using  it  to  generate  synthetic  data. 
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CHAPTER  IV: 

SYNTHETIC  DATA  GENERATION 


IV.l)  Overview 

The  generation  of  high  quality  synthetic  pressure  fields  is  an  important  branch  of  acoustic 
research.  Because  present  methods  can  only  approximately  compute  the  fields  associated  with  a 
point  source  for  complex  environments,  simplified  environments  are  often  considered  for  which 
fewer  approximations  must  be  made.  One  important  environment  which  lends  itself  well  to 
analysts  but  which  has  sufficient  complexity  to  be  of  interest  for  practical  problems,  is  the  horizon¬ 
tally  stratified  environment.  It  is  an  excellent  model  for  the  conditions  present  in  the  deep  ocean 
over  an  abyssal  plain,  and  consequently  of  direct  interest  to  us.  Currently,  techniques  for  calcu¬ 
lating  synthetic  fields  arising  from  a  CW  point  source  in  this  environment  exist  in  the  literature 
[7, 11,3  ]  These  techniques  are  based  upon  the  numerical  evaluation  of  the  Sommerfeld  integral  [1 
J,  for  which  two  major  computational  efforts  are  required.  First,  the  plane  wave  reflection 
coefficient  for  the  bottom  profile  must  be  numerically  generated.  For  this  the  propagator  matrix 
method  [12,8  ]  is  used.  Second,  the  pressure  field  is  computed  as  the  Hankel  transform  of  the 
depth-dependent  Green's  function  (which  is  simply  derived  from  the  plane  wave  reflection 
coefficient).  Typically,  in  these  techniques  many  of  the  degradations  associated  with  the  numeri¬ 
cal  evaluation  of  the  Hankel  transform  are  not  carefully  controlled.  In  this  chapter  we  exploit  the 
properties  of  the  Hankel  transform  derived  in  Chapter  II  to  carefully  control  these  errors.  We 
will  show  in  Section  (IV  3a)  that  a  major  source  of  numerical  error  is  aliasing,  which  becomes 

important  because  asymptotically  the  fields  decay  only  isy  .'  We  associate  this  slow  rate  of  deca-, 

with  the  source  singularity,  — -  In  the  depth-dependent  Green’s  function  and  show  how 

WHk? 

to  separate  this  portion  of  the  computation  from  the  numerically  computed  Hankel  transform. 

1)  Asnmiai  tor  the  present  that  there  ue  no  tripped  modes  ueocieted  with  low  speed  layers  within  the 
bottom. 
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The  remaining  numerical  calculations  are  significantly  less  degraded  by  aliasing  and  are  well 
behaved  in  general.  They  remain  in  the  form  of  a  Hankel  transform  for  which  we  can  exploit  the 
computation  efficiencies  now  available  [9  ]  The  result  is  a  new  hybrid  procedure  which  is  compu¬ 
tationally  efficient  and  significantly  more  accurate  than  existing  methods.  This  hybrid  scheme  is 
illustrated  with  examples  of  synthetically  generated  fields. 

In  Section  (IV Jb)  we  discuss  the  difficulties  associated  with  numerically  evaluating  the 
Hankel  transform  of  the  depth-dependent  Green’s  function  when  stow  speed  layers  are  present  in 

the  bottom  which  give  rise  to  proper  modes.  Proper  modes  are  associated  with  the  — — - 

singularities  (with  kr  near  the  real  axis)  in  the  depth-dependent  Green’s  function  and  contribute 

terms  to  the  field  which  decay  asymptotically  as  — \s .  This  very  slow  decay  in  the  field  causes 

Vr 

severe  aliasing  problems  when  it  is  calculated  using  a  numerical  Hankel  transform  algorithm.  In 
Section  (IV .3b)  we  show  how  to  separate  the  effects  of  proper  modes  from  the  numerical  calcula¬ 
tions  by  performing  part  of  the  Hankel  transform  analytically.  We  make  this  separation  in  such  a 
manner  that  the  portion  of  the  field  assigned  to  the  analytical  calculations  is  exact  and  finite  for 
ali  ranges.  This  makes  it  possible  to  numerically  calculate  the  residual  numerical  contribution  to 
the  field  accurately  and  add  the  result  to  the  analytically  determined  expression.  The  result  of  the 
total  hybrid  algorithm  is  a  field  which  is  accurate  for  all  ranges  and  which  can  accurately  include 
the  effects  due  to  proper  modes  arising  in  the  presence  of  slow  speed  layers.  We  present  an  exam¬ 
ple  of  a  field  generated  synthetically  with  this  total  hybrid  method  and  show  how  the  result  is 
significantly  better  than  what  would  have  been  achieved  without  removing  the  effect  of  the  poles. 

In  this  chapter  we  also  develop  a  numerical  implementation  of  the  propagator  matrix  method 
for  generating  the  plane  wave  reflection  coefficient  that  is  well  behaved  numerically.  We  begin 
the  chapter  by  describing  the  computation  of  the  plane  wave  reflection  coefficient  by  means  of  the 
Thomson-Haskell  method  [12,8  ]. 


F 
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IV.2)  The  Propagator  Matrix  Approach  to  Generating  the  Plane  Wave  Reflection  Coefficient 


a)  The  Method  in  Principle 


I)  Overview 

To  calculate  the  plane  wave  reflection  coefficient  we  consider  the  response  of  a  layered  bot¬ 
tom  to  an  incident  plane  wave  as  shown  in  Figure  lV.2a.i.l.  Within  the  n'h  isovelocity  layer  we 
express  the  field  as  the  vector: 


/» (*) 

l/(*)j/o(*'r)e  (1) 

where  PH(z)  is  the  pressure  in  the  n‘h  layer  and  UK{z )  is  the  normal  component  of  the  velocity. 
We  have  chosen  this  representation  because  P(z)  and  U (z )  are  continuous  in  z ,  even  across 
layer  interfaces.  In  the  discussion  which  follows  we  will  suppress  the  time  and  radial  dependence 
of  the  field  because  they  are  the  same  in  all  layers. 

In  the  propagator  matrix  approach,  the  impedance  at  the  bottom  layer: 

r  -  )  r) 

is  available  from  the  material  parameters.  In  principle  this  impedance  is  used  to  determine  the 
reflection  coefficient  at  the  top  interface  in  three  steps.  First  the  field  at  the  top  interface  is 
related  to  the  field  at  the  bottom  interface  by  the  propagator  matrix: 

P(z  o)j  [*(*»)  I 

l/(z  o)|  =  <t>  [l/(rv)j  (3) 

Next  the  incident  and  reflected  pressure  waves  at  the  surface  are  related  to  the  field  at  the  top 


interface  and  then  to  those  at  the  bottom  through: 

>  *.o]  [/*(*<>)]  [?(**)]  [«11  «2l]  p*(*w) 

P  -,0  U  (Zy)  (Xu  «22  U  (*Ar) 

Finally,  the  reflection  coefficient  is  calculated  in  terms  of  the  impedance,  ly+i,  using 


P+,0  ®U  *21  1 

P -,0  *21  *22  CjV+l 


i 
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Figure  lV.2a.;.l  The  waves  generated  in  a  layered  bottom  in  responses  to  an  incident  plane 
wave 
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so  that 


P-.0  _  Oil  +  £n  -*-l“l2 
P+fl  “21  +  £tf+l“22 


(6) 


11)  The  Propagator  Matrix 

The  essential  element  of  this  approach  is  the  propagator  matrix,  <t>,  of  Equation  IV.2a.13.  In 
this  section  we  review  its  derivation. 


Within  any  isovelocity  layer,  the  field  can  be  considered  as  the  superposition  of  a  positive 
and  a  negative  traveling  wave.  The  pressure  field  is  given  by 

P(*)»P+e*’‘  +  P-e~*'’  (1) 

The  normal  component  of  velocity,  U  (2  ),  is  related  to  P  (r  )  through  the  telegraph  equations  [2  ] 

For  the  non-normal  case  we  use  — - —  “  m  —itapU  which  implies  that: 

dt  dt 


U(z)  •  r —p+e'*’*  - 

1  up  tup 


or  defining  Y  0 


up 


1/(2)  -  YoP+e*1*  -  YoP-e 
In  matrix  form  Equations  (1)  and  (3)  become: 


~ik,i 


P(l) 

U(t) 


ik,l 

Y* 


e  ’  e 

ii,l  „  ~ik,. 


If 


X*l) 

u(*0 


is  known  at  some  point  in  the  layer  then 


-ik,l 

Y* 

X*  2) 


P  + 

P  _ 


(2) 


(3) 


(4) 


inverting  Equation  (4)  to  find 


P  + 
P_ 


in  terms  of 


X*  1) 


can  be  computed  in  principle  by 

and  then  calculating 


Pit  2) 

IX*2>. 


from 


P* 

P. 


Combining  these  operations  into  one  step  gives: 


X*  2) 

1/(22) 


-iti, 

e  •  •  e 


Y* 


t*,i, 

Y* 


e  ■"*  e 

^1*1  _y  - 


-I'd*' 


-1 

>(*l) 

V  (2t) 

(5) 


which  gives  us: 
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PM 

I'M 


cos  *,(2  2-2  0  -rr-sinfc»(*2-:i) 

*  0 


X*  l) 

U(z  i) 


(6) 


jiYosin^kj-Jt)  cos^,(r2-zt) 

The  values  of  k,  and  Y  0  are  functions  of  the  material  parameters  of  the  layer  under  considera¬ 
tion.  In  particular  if  cn  is  the  sound  speed  in  tayer  n,  p„  its  density,  kr  the  horizontal  wave 
number  of  the  incident  plane  wave  (by  Snell’s  law  common  to  all  layers),  and  to  the  temporal  fre¬ 
quency  of  the  CW  source,  then  kH  ■  — ,  kt  ■  \/k}-k?,  and  K0  » 

c„  tap 

To  indicate  explicitly  the  dependence  on  the  material  properties  of  the  nlh  layer  we  write: 


X*  2) 

V  (zj) 

when  2 2  and  zt  are  both  within  layer  n 


X*i)] 

V{*  i) 


(7) 


To  calculate  the  field  at  the  top  interface  in  terms  of  the  field  at  the  bottom  interface,  as 
shown  in  Figure  IV2ai.l,  we  can  use  the  previous  discussion  which  was  applicable  only  to  a  single 
layer,  to  relate  the  field  at  z„  to  the  field  at  z*.^: 


>(*„-<) 

UMx) 


♦n(z._!-z,) 


XO 

UM 


We  then  iterate  the  procedure  through  ail  the  layers  to  find 


X*  o) 
V{z0) 

b)  Numerical  Implementation 


\*l 


P(Zn) 

UM 


PM 

UM 


PM 

VM 


(8) 


(9) 


I)  The  modified  propagation  matrix 

The  bulk  of  computation  associated  with  the  propagator  matrix  approach  is  the  accumulation 
of  the  matrices  When  these  are  accumulated  on  the  computer,  the  actual  opera¬ 

tion  is  ♦i(<l>2(  '  '  ’  $#)))•  ^  >*  pottiMe  for  the  scale  of  the  accumulated  product  to  differ 
dramatically  from  any  particular  <t>( .  Because  of  the  limited  dynamic  range  in  the  computer  it  is 
advisable  to  scale  terms  to  make  them  comparable  before  accumulation.  Fortunately  the  final  cal- 
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culation  for  the  reflection  coefficient  depends  only  on  ratios  of  elements  in  4>.  For  this  reason  we 
normalize  each  of  the  so  that  its  largest  value  equals  1.  This  procedure  alone  could  cause 
another  problem  stemming  from  the  different  scales  in  general  for  P  and  U ,  which  is  due  to  their 
different  units.  To  bring  P  and  U  into  the  same  units  we  do  not  actually  relate 


>(*.-1) 


to 


>(*.) 

UM 


(i) 


but  rather  consolidate  units  by  multiplying  the  normal  component  of  velocity  by  the  characteristic 
impedance  of  that  layer.  Therefore  we  actually  calculate: 


where 


>(*o) 

a\  *1 

aV  bN 

5otf(*o), 

m 

tib\  £iai 

knbN  htaN 

iNuM 

"  vU/c, 


at  m  co*  4  (*(  ~zi  -1) 

bt  ■  -isin  1) 


U)  Relation  of  the  modified  propagation  matrix  to  the  incident  and  reflected  waves 


(2) 


We  now  relate  the  field  variables  to  the  incident  plane  wave  and  the  resulting  reflected  plane 
wave  by  slightly  modifying  Equation  (lV^aii.4).  We  assume  that  the  top  interface  is  at  2  -  0  so 
that: 


and 


P(  0) 

1  0 

5o  1/(0) 

0  to. 

Fo 


1  1  ' 

P*. 0 

Y0  -Yo 

P-fi 

1  0 

P(0) 

0  yo 

ioU(0) 

(1) 


(2) 
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By  defining: 


4>u  4»ia 

N 

at  bt 

*21  4*22.. 

-  n 

i-l 

hibi  Sifl, 

(3) 


and  using  Equation  (IV.2b.i.2)  we  have: 


P+fi 

1 

1  1 

4>n  4>i2 

PM 

P-fl. 

as  — 

2 

1  -1 

4>21  4>22. 

S  NU(zs) 

(4) 

We  now  need  to  use  the  fact  that  the  pressure  and  velocity  fields  in  the  last  layer  are  made  up  of 
only  positive  traveling  waves  so  that  (referring  to  Equation  IV.2aii.4)  PN+l**P(zN)  and 
1 


UN, 


i 


P  (?N )  we  have 


If  we  now  use 


>♦.  0 

m± 

4>tl+<h21  4*52  "^4*22 

/-.o 

2 

4>u_4>2i  4>u_4>22 

£w+l 


N  +1 


(5) 


+1 


Jy_ 


we  have  the  reflection  coefficient 


(6) 


„  P  -,o  4>n-«h2i  +  +i(<hi2"‘J>2z) 

r  ■  ■  ■’  m . .  — i ■  —i  . -  (7) 

P+.o  <hn+$2i  +  --i(4>i2+4>22) 

Equations  (IV-2bii.2)  and  (7)  show  that  this  approach  uses  only  the  ratios  of  the  impedances 
in  adjacent  layers  and  never  the  impedances  themselves.  These  ratios  are  much  better  behaved  in 
general  than  the  individual  impedances.  For  this  reason,  because  the  use  of  P  and  YqU  instead  of 
P  and  U ,  and  because  of  the  scaling  of  the  layer  propagation  matrices  this  implementation  of  the 
propagator  matrix  approach  has  good  numerical  properties. 


c)  Selected  Properties  of  the  Reflection  Coefficient 

In  Figure  IV.2C.2  we  present  a  perspective  plot  of  the  log  magnitude  of  the  reflection 
coefficient  as  a  function  of  kr  tor  the  bottom  of  Figure  IV ,2c.  1  calculated  using  the  numerical 
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f  =  220  Hz 
|z  +  z0  |  =  2m 
k0=. 8975979  m-1 


C0  s  1540  m/s 
po  -  1  g  /cm3 

C,  *  1493.8  m/s 
yo,  *  1.5  g/cm3 

C2  *  1700  m/s 
pz  =  2.0  g/cm3 


Figure  IV.2c.l  The  bottom  parameters  used  to  generate  the  reflection  coefficient  shown  in  Fig¬ 
ure  lV.2c.2 


algorithm  just  described.  The  reflection  coefficient  is  displayed  both  for  those  horizontal  wave 
numbers  corresponding  to  real  angles  of  incidence  (0  <  Rt(kr )  <  k0  where  k0  is  the  water  wave 
number)  and  for  horizontal  wave  numbers  corresponding  to  complex  angles  of  incidence 
(k  o  <  Rt(kr )).  The  complex  angles  correspond  to  evanescent  waves.  Single  evanescent  waves  do 
not  carry  any  time  average  power  flow  (their  Poynting  vector  is  imaginary)  and  consequently  the 
magnitude  of  the  reflection  coefficient  is  not  limited  to  be  less  than  one  in  the  evanescent  region, 
k0<kr  [6] 

In  Figure  1V.2c2  a  pole  is  apparent  in  the  reflection  coefficient  on  the  real  kr  axis  in  the 
evanescent  Tegion.  This  on  axis  pole  corresponds  to  a  proper  mode  propagating  in  the  low  speed 
layer  within  the  bottom.  Other  off  axis  poles  corresponding  to  leaky  modes  are  also  apparent  in 
the  reflection  coefficient.  A  discontinuity,  or  cut,  can  be  seen  extending  from  k2  along  the  real  kr 
axis  to  infinity.  This  is  the  branch  cut  extending  from  the  branch  point  at  Jfcj.  Another  cut  extend¬ 
ing  from  k0  to  infinity  falls  on  this  same  line  and  is  therefore  not  apparent.  The  origin  of  these 
branch  points  and  cuts  can  be  found  in  our  derivation  of  the  reflection  coefficient  where  we  asso¬ 
ciated 

eiVki-k*i  with  p +  (1) 

and 

e-iVkPi ?.  with  p  _  (2) 

Clearly  the  roles  of  P  +.  and  P  _  would  be  reversed  by  changing  the  choice  of  sign  for  the  square 
root.  For  incident  and  reflected  wave  this  would  correspond  to  inverting  the  reflection  coefficient 
(if  no  other  waves  were  affected).  The  two  sheets  corresponding  the  the  branch  point  at  /fc0 
reflect  the  two  choices  of  sign  for  the  incident  wave.  We  have  displayed  the  choice  associated  with 
positive  real  power  flow  for  the  incident  wave. 

In  the  intermediate  layers  such  as  layer  1  of  this  example,  changing  the  role  of  P  +  and  P  _ 
would  not  affect  the  reflection  coefficient.  For  the  intermediate  layers,  the  physical  problem  does 
not  name  (or  distinguish  between)  forward  and  backward  traveling  waves.  Consequently  there  are 


'  s-**  !'.«*** 
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no  branch  points  associated  with  intermediate  layers. 

t  n  •  1  q 

If  the  opposite  sign  were  chosen  for  the  square  root,  V^v+i  Hfc,*,  associated  with  the 
emerging  wave,  F+jy+i  front  the  bottom  of  the  stack  of  layers  (into  the  isovelocity  half  space)  the 
direction  of  energy  flow  associated  with  that  wave  would  change.  Unlike  the  intermediate  layer, 
there  is  no  returning  wave  in  the  isovelocity  half  space.  Consequently,  the  physical  problem  would 
change.  For  this  reason  we  see  a  branch  point  at  kv+ 1  reflecting  the  two  different  "physical''  prob¬ 
lems. 

In  Figure  IV.2c.2  we  have  chosen  to  display  the  Riemann  sheet  for  which  both 
Re(Vko  ~*7)  >  0  and  Re(\4jv+i  “*/)  >  0-  On  this  sheet  only  waves  with  real  power  flow  in 
the  positive  direction  are  associated  with  P  +.  This  constrains  our  incident  waves  to  be  those  with 
power  flow  from  the  source  to  the  layered  oottom  and  specifies  that  there  is  no  power  flow  return¬ 
ing  from  infinity. 

When  we  perform  the  integrations  discussed  later  we  must  choose  which  side  of  the  cuts  to 
integrate  upon.  For  reasons  of  convergence  we  choose  the  side  for  which  Im(\A/~V)  >  0 
when  j  =  0  and  N  +1,  is  satisfied.  Consequently,  whenever  we  integrate  the  reflection 
coefficient  in  the  complex  kr  plane,  we  always  satisfy  both  >  0  and 

Im(‘ S/k  J-k?)  >  0  for  ;  -  0  and  N  +1. 

IV .3)  Evaluating  ike  Sommerfeld  Integral 

Once  the  plane  wave  reflection  coefficient,  T (Jkr ),  has  been  computed  it  is  necessary  to 
evaluate  the  Sommerfeld  integral: 


P*(r)  -  /  — r(fc, )e ‘ 0(rkr )kr dkr 

n  v*- 1 


(1) 


in  order  to  compute  the  reflected  pressure  field.  The  Sommerfeld  integral  is  in  the  form  of  a 
Hankel  transform  of  the  depth  dependent  Green’s  function. 


G(krjj  0) 


v£F*7 


IX*,)* 


i»  ■rial 


The  properties  of  the  Hankel  transform  were 


I 


developed  in  Chapter  11.  In  Sections  (6)  and  (7)  of  that  chapter  we  discussed  the  effect  of  trun¬ 
cating  the  integration  at  some  finite  value,  and  the  effect  of  sampling.  The  truncation  was  accom¬ 
plished  by  multiplying  the  function  to  be  transformed  by  some  finite  length  window.  For  the  gen¬ 
eration  of  synthetic  data  we  find  that  windowing  of  the  Green’s  function  is  not  an  important  con¬ 
sideration  in  general  because  when  z  +z0  >  0,  G  ( kr,z  ,:0)  decays  exponentially  in  kr  for  fc,  >  kQ. 
Except  when  z  +io  is  very  small  we  can  integrate  Equation  (T)  until  the  Greens  function  is  negli¬ 
gible  and  truncate  at  that  point.  It  is  not  necessary  to  multiply  by  a  windowing  function. 

The  issue  of  sampling  and  the  associated  degradation  introduced  into  the  transform,  aliasing, 
can  be  very  much  a  problem  however. 


a)  The  Source  Singularity 

In  order  to  highlight  the  issues  associated  with  the  source  singularity, 


r,  and  the 


vSF*?' 

propagation  terms  in  the  Sommerfeld  integral,  we  first  consider  the  evaluation  of  Equation 
(IV 3.1)  for  a  hard  bottom  case  where  T (fcr )  =  1. 

For  the  hard  bottom  case  the  pressure  field  is  given  by  the  known  integral: 


i\ZkfT2it+  i.i 


Joilrftdt 


ii0V»-2+(»+»o)2 


(2) 


o  .  v^+cT+731 

We  evaluated  this  integral  numerically  with  the  Fourier-Bessel  series  to  obtain  an  estimate  for  the 
field: 


7  N  ®  (~T,Z 

'*”■7? - -  0<'<A  (3) 

where  for  this  example1,  A  was  chosen  to  be  2000  and  I  *  +ro 1  “  2  .  In  Figure  (IV3a.l)  we  com¬ 
pare  the  log  magnitude  of  the  result  (dots)  with  the  known  transform  (solid  curve).1  We  see  that 
the  magnitude  of  the  numerically  generated  held,  PR(r ),  oscillates  rapidly  in  contrast  with  the 


l)The  output  of  the  Fourier-Bessel  series  be*  been  displayed  to  twice  its  region  of  validity  to  better  illus¬ 
trate  the  source  of  degradation. 


i 
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ue  Hankel  transform.  As  we  will  now  show  these  oscillations  are  due  to  aliasing  in  the  numeri¬ 
cally  computed  Hankel  transform. 

In  Section  (11.7)  we  showed  that  the  effect  of  sampling  on  the  Hankel  transform  is  to 
approximately  produce  an  aliased  estimate  of  the  true  transform,  VrPs(r).  Since  for  this  exam¬ 
ple,  Pr  ( r )  decays  asymptotically  as  1/r ,  ^tPr  (r )  decays  asymptotically  only  as  1/Vr* .  What  we 
see  in  Figure  (IV.3a.l)  is  given  approximately  by: 

ik^+b+ttf  _  a0V(2a-r^(,>,0)J 

0<  r  <  2A  i>»(r)l  «  -p  - -  =r  -  - = . =  =  (4) 

^  VrT+?T+^j!  V(2A  -r  )z+(*  +*0)2 


When  r  is  much  greater  than  z  +z<>,  (r )  is  approximately: 


1-  ii </  Ucal2A-r  I  1  I 

-  ^r\\ 


Since  we  are  in  the  region  r  <  2A  this  can  be  rewritten: 


11  0t(/  UtoH  ]| 

II 

We  can  write  Equation  (6)  in  terms  of  the  desired  transform  and  a  modulation  term  as: 

1  ff  ik0 2A  uk„2A 

—=  |  —p — t. — = —  e  °r  +  2t- -aagagsin  kf/ 


Which  upon  defining  «(r  ) 


!/>*(*) I  «  e‘V  +2‘'^«n*0rj  (8) 

When  r  «  2A  e(r)  is  small,  so  that  the  magnitude  of  PR  (r )  appears  as  roughly  the  correct 
transform  with  a  modulation  term. 


We  note  at  this  point  that  if  we  had  sampled  the  output  of  our  transform  at  an  integral  mul¬ 
tiple  of  2rr/ko  we  would  not  have  seen  these  oscillations.  At  this  sampling  rate  the  cosine  would 
have  appeared  as  a  DC  offset  in  the  magnitude  of  the  pressure  field.  If  the  output  sampling  rate 
were  near  but  not  exactly  an  integral  multiple  of  2rr/k0  the  cos(k( /)  would  have  appeared  as  a 
low  frequency  modulation  because  the  sampling  is  in  effect  demodulating  the  cosine  down  to  a 
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low  (but  not  now  zero)  frequency.  This  result  is  an  important  one  because  frequently  pressure 
fields  are  generated  by  using  an  FFT  based  approximation  to  the  Hankei  transform  (IV.3.1)  and 
the  water  wave  number  is  the  maximum  wave  number  used  [5,4  ]  The  grid  resulting  from  such 
processing  is  exactly  an  integral  multiple  of  2ir  /k0.  Carrying  the  integration  to  higher  wave 
number  would  make  evident  the  modulation  in  the  answer  by  automatically  providing  the  output 
on  a  finer  grid. 

The  problem  of  aliasing  arose  because  the  field  being  computed  decayed  only  as  —  which 

r 

forces  us  to  use  a  very  high  sampling  rate  to  properly  sample  the  Hankei  transform.  We  now  note 
that  this  —=  decay  is  due  to  the  — . . —  singularity  in  the  Green’s  function.  It  is  well  known 

Vr  WHk? 

that  the  asymptotic,  or  far  field,  character  of  such  a  transform  is  determined  by  the  singularities 
of  the  kernel  over  the  path  of  integration  [10  ].  The  Green’s  function  which  is  transformed  in 
Equation  (IVJ.l)  was 

j- /».  _  _  \  _  i  rn. 

G<‘’”0)  <9) 

The  asymptotic  character  of  the  transform,  P  (r ),  is  dominated  by  the  singularity 


/k.Q  ~kr2 


The  integral 


77 =F=r  -  /  j—  S'fiFZ" ' J0(k,r)krdkr  (11) 

VrV  o  y/k o  -kr 

shows  us  that  this  singularity  is  in  fact  associated  with  the  1/r  decay  rate.  Physically  this  singular* 
ity  was  due  to  the  angular  spectrum  of  the  point  source.  The  1/r  decay  associated  with  this  singu* 
larity  is  often  associated  with  the  point  source  by  noting  that  the  field  around  a  point  source  must 
decay  at  that  rate  in  a  manner  such  that  the  intensity,  which  decays  as  the  field  squared, 
integrated  over  any  three  dimensional  shell  enclosing  the  point  source,  would  not  be  a  function  of 


r. 


The  source  of  stow  asymptotic  decay  we  have  isolated  suggests  a  procedure  for  reducing  the 
problem  of  aliasing.  We  remove  the  source  singularity  from  the  kernel,  numerically  transform 
what  is  left,  and  add  the  result  to  the  analytically  determined  transform  of  the  source  singularity. 
When  we  remove  the  singularity  we  must  do  so  in  a  manner  such  that  the  numerical  transform  we 
must  perform  is  well  behaved.  We  allow  for  a  general  T (Jkr )  but  at  this  time  assume  that  T(kr) 
has  no  singularities  along  the  real  kr  axis  with  asymptotic  contributions  to  the  field  capable  of 


dominating  those  of  the 


i 

VkiFi? 


singularity. 


To  this  end  we  rewrite  integral  (IV .3.1)  as: 


PK(r)~S 


i _  i\£F£?n+«oi 


Jo(krr)krdkf 


'/  [r(fc,)-r(fc )]  ^  t  V^,,4,o'j0(jt,r  )KdK  +  T(*)jf  ^=p«,V^',+,°Vo(*,r)*,dtr  (13) 


If  we  define: 


L{kr)  m  [r(*r)  -  r  (k)]--^— (14) 

Vk0  ~k, 

so  that  L  (kr )  does  not  have  the  1  /'S/k  J —k}  singularity  at  kr  =  k  q.1  then  we  can  write  Equation 


(13)  as: 


ii  V>  *+(»+,„ 


P»(r)  =  /  L(kr)J0(krr)krdkr  +  r(*)~===  (15) 

o  V>2+(r  +r0)2 

Because  L{kr)  does  not  have  this  singularity  along  the  path  of  integration  the  output  of  the 
numerical  transform  will  decay  at  a  rate  faster  than  1/r.  The  asymptotic  1 It  decay  is  provided  by 


We  show  in  the  appendix  that  if  the  impedance  and  its  first  derivative  at  the  interface  is  finite  for 
kr  =  ko  then  the 


lim  L(kr ) 
*o 


Z  i 
“Po 


where  Z\  is  the  impedance  of  the  bottom  at  *,  «  k o,  u>  is  2ir  source  frequency,  and  po  is  the  density  of 
the  water.  For  an  isovelocity  half  space  this  expression  reduces  to 


L{k  o) 


Pi 

povTFiJ 


Which  is  finite. 


the  analytic  term  which  can  be  recognized  as  the  specular  reflection  when  r  is  very  large  (glancing 
incidence).  These  observations  are  confirmed  in  the  examples  which  follow. 

In  the  following  examples  we  illustrate  the  generation  of  synthetic  pressure  fields  through 
the  hybrid  algorithm  implied  by  Equation  (15)  where  the  integral  is  performed  with  a  numerical 
Hankel  transform  algorithm  and  the  analytical  expression  is  the  result  of  integrating  the  singular¬ 
ity. 

I)  Hard  Bottom 

This  is  the  degenerate  example  because  for  T(kr )  constant,  the  entire  transform  is  per¬ 
formed  analytically.  The  result  of  the  analytic  transform  was  compared  to  the  direct  numerical 
transform  Figure  (IVJa.l). 

II)  Slow  bottom 

Figure  (IV.3aai.l)  shows  the  bottom  parameters  for  this  example.  Figure  (IV3aii.2)  shows 
the  result  of  the  hybrid  calculation  (solid  line)  versus  a  direct  numerical  calculation.  The  improve¬ 
ment  is  dramatic.  Figure  (IV.3aii3)  compares  the  hybrid  field  of  Figure  (lV.3aii2),  with  its 
numerically  generated  component.  As  can  be  seen,  the  near  field  is  dominated  by  the  numerically 
generated  component.  As  range  increases  this  numeric  term  begins  to  suffer  from  aliasing  prob¬ 
lems  but  the  analytic  term  begins  to  dominate,  minimizing  the  effect  of  aliasing  on  the  computed 
field  at  large  ranges. 

ill)  Fou  Bottom 

Figure  (lV3aaii.l)  shows  the  parameters  of  the  fast  bottom  for  this  example.  Figure 
(IV3ajii2)  shows  the  hybrid  calculation  versus  the  direct  numerical  calculation.  Figure 
(IV3aiii3)  presents  the  hybrid  field  and  its  numeric  component.  The  improvements  are  similar  to 
the  fast  bottom  case. 


b)  Poles  Due  to  Slow  Speed  Layers 
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f  =  220  Hz 
|  z  +  z0|  =  2m 

k0=. 8975979  m“’ 


Co  s  1540  m/s 
p0  a  I  g/cm3 

'7  >/ '/////■ -7/  ■////./ 

C,  =  1493.8  m/s 
y0,  s  1.5  g  /cm3 


Figure  TV.3a.il.l  Parameters  of  bottom  used  for  slow  bottom  example 
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Figure  IV.3«.H.3  Comparison  of  hybrid  field  calculated  for  slow  bottom  example  and  its  numeri 
cally  generated  component 
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C0  s  1 540  m/s 
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Ci  *1700m/s 
P\  s  2.0  g/cm3 
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Figure  (IV 3b. 1)  shows  the  parameters  of  a  slow  speed  layer  between  two  isovelocity  half 
spaces.  These  are  the  same  parameters  used  to  generate  the  perspective  plot  of  the  reflection 
coefficient  presented  in  Figure  (IV .2c .2).  Figures  (IV.3b.2a)  and  (IV 3b 2b)  show  the  magnitude 
and  phase  of  the  reflection  coefficient  for  this  bottom  as  a  function  of  horizontal  wave  number. 
We  see  that  for  this  example  the  reflection  coefficient  has  a  singular  point  beyond  the  water  wave 
number,  k0.  That  singularity  is  a  simple  pole  associated  with  a  proper  mode  excited  in  the  low 
speed  layer.  Such  a  proper  mode  can  appear  only  for  fco  <  kr  <  fcy+1.  In  this  region  conserva¬ 
tion  of  energy  is  not  violated  because  the  waves  are  evanescent.  Proper  modes  are  generated 
when  the  low  speed  layer  acts  like  a  dielectric  waveguide.  When  this  happens  energy  diffuses  (tun¬ 
nels)  from  the  point  source  to  the  low  speed  layer  but  does  not  otherwise  propagate  vertically. 
£nergy  from  the  field  is  now  constrained  to  decay  in  only  two  dimensions  instead  of  three  and  we 

expect  that  the  field  associated  with  the  pole  will  decay  asymptotically  as  — so  that  the  integral 

vr 

of  the  flux  over  any  two-dimensional  ring  surrounding  the  source  remains  constant. 


Poles  such  as  this  disrupt  the  asymptotic  character  of  the  field  derived  in  the  previous  sec¬ 
tion.  As  before  we  would  like  to  analytically  determine  the  contribution  of  these  poles  and 

i 


remove  them  as  we  removed  the 


integral: 


t(rj+t0tkr)  m  J 


V*o 


singularity.  To  do  so  it  is  necessary  to  evaluate  the 


*2-*  2 
*fTi 


i  +i0 


Vfc2-*,1 


J0(krr)krdk, 


(1) 


In  Appendix  (I)  we  show  that  for  Irti(k,()  0  (associated  with  no  return  from  r  *  *>  ) 
l(r  j+z0\kr  )  is  given  by: 


II+I01 


(2) 


where 


0  -  +\/k?-ki 


(3) 


i 


f  =  220  Hz 
|  z  +  z0  |  =  2  m 
k0*. 8975979  rtT1 


C0  =  1540  m/s 
p0  =  1  g  /cm3 


C,  =  1493.8  m/s  * 
yo,  =  1.5  g/cm3  I 


10  m 


C2  =  1700  m/s 
p2  -  2.0  g/cm3 


Figure  IV.3b.l  Parameters  of  bottom  used  for  slow  speed  layer  example 
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The  first  integral  is  easily  evaluated  on  the  computer.  In  addition  as  r  becomes  large  the  first 
term  rapidly  approaches: 


i‘±oV>1+(i+io)J 

The  second  term  decays  as  l/v7  when  Jkf(  is  real. 


(4) 


Equation  (2)  is  also  correct  for  Im(k,)  >  0,  but  when  Im(fc, )  »  0  the  poles  no  longer 
contribute  asymptotically  as  1  /V7  because  the  Hankel  function  decays  exponentially.  Under  these 
conditions 


^  (*r, )  1*,  ■' 

< — =— e  e  ‘ 

v7 


(5) 


(6) 


As  Im(jfcr )  becomes  large  the  exponential  decay  dominates  the  1/vT  decay  even  over  the  finite 
range  that  concerns  us.  It  is  for  this  reason  that  we  consider  only  those  poles  near  the  real  axis 
(close  to  the  path  of  integration)  and  leave  the  others  to  the  numerical  part  of  the  transform. 

With  /  (r  ,z  +Zq^)  *o  defined,  the  reflected  pressure  field  can  be  written: 


MO-J 


Vkif-k} 


r(*')-27drr 


j^ktf)kfdkr  +  ^^/(r  ,*  +Io;  *,)  (7) 


Where  the  expression  in  brackets  no  longer  has  any  poles  near  the  line  of  integration  and  so  can 
be  evaluated  as  before. 

In  order  to  remove  the  poles  as  required  in  Equation  (7)  it  is  necessary  to  determine  with 
precision  the  pole  locations,  kr.,  and  their  scales,  (a_;)( .  The  pole  locations  can  be  found  using 
standard  complex  root  finding  techniques,  though  care  must  be  taken  to  provide  the  root  finding 
algorithm  with  values  of  the  reflection  coefficient  on  the  Riemman  surface  so  that  it  appears  ana¬ 
lytic  except  at  isolated  singularities.  This  means  that  the  branches  chosen  for  the  square  roots 
must  be  taken  in  such  a  manner  that  a  branch  cut  is  never  placed  between  points  simultaneously 
considered  by  the  root  finder.  Once  the  root  locations  are  known,  the  scale  factors  can  be  found 


for  those  singularities  far  from  any  others  by  determining  a  least  squares  fit  to: 


r<V  -  7 T7i 


j  =  1,2. 


,N 


(8) 


provided  that  the  kr^  are  taken  sufficiently  close  to  kTi  that  T(kr )  is  well  approximated  by  just  one 
pole  in  that  region. 

If  many  poles  are  clustered  together,  they  can  be  determined  simultaneously  by  solving: 


r(V=  .2 


(a-t)t 


j  =  1,2, 


,N 


(9) 


kr^-kr,2 

for  N  sufficiently  large.  If  a  pole  is  near  a  branch  cut  then  the  poles  on  the  other  side  of  the  cut, 
on  the  opposite  sheet,  and  near  the  cut  must  also  be  considered  to  be  near  that  pole. 


Figures  (IV ,3b 3a)  and  (IV3b3b)  show  the  magnitude  and  phase  of  the  reflection  coefficient 
of  Figure  (IV 3b .2a)  minus  the  pole  contribution: 


Fot 


this 


example 


!•(*,)- 


a-l 

(*r2-*r.2) 


a_!  -  1.689712*  10"1 


kr  m  9.069830*  10_1  +  i  2.488749*  10~3. 


5.027826*  10-4 


(10) 

and 


Note  the  difference  in  scale  between  Figures  (IV3b2a)  and  (IV 3b 3a).  The  small  notch  visi¬ 
ble  at  kr  ~  kT[  is  due  to  a  small  amount  of  error  in  the  estimate  of  kr< . 

A  notable  feature  of  Figure  (IV  ,3b  .3a)  is  the  unmasking  of  the  off  axis  zeros  in  the  region 
k/t+i  <  kr  <  k  o  where  previously  I T  (k, )  I  =■  1.  These  zeros  can  be  clearly  seen  in  the  perspec¬ 
tive  plot  of  the  total  reflection  coefficient  in  the  complex  plane  that  was  presented  in  Figure 
(1V2c2). 

Figure  (IV3b.4)  presents  the  hybrid  field  (solid  line)  versus  the  field  calculated  without 
removing  the  pole  from  the  reflection  coefficient  (but  otherwise  removing  the  l/\/kJ-k?  singu¬ 
larity  as  in  the  previous  hybrid  examples).  The  spread  in  the  directly  computed  field  due  to  alias¬ 
ing  is  severe  because  aliasing  in  the  Hanksl  transform  severely  affects  a  function  that  decays  as 


MAGNITUDE 
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1/Vr  .  The  hybrid  field  does  not  exactly  follow  the  contour  of  the  top  of  the  spread  just  as  the 
hybrid  computations  in  the  previous  examples  did  not  exactly  follow  those  contours  when  (he 
aliasing  became  severe.  Figure  (IVJbJ)  presents  the  log-magnitude  of  the  analytically  generated 
pole  contribution  (solid  line)  and  the  remainder  of  the  field  exclusive  of  the  pole  contribution. 
The  non-pole  contribution  is  most  significant  for  short  ranges,  while  for  this  near  bottom 
geometry  the  pole  contribution  dominates  farther  out. 

The  expression  for  /  ( r  ,z  +Zo’,kr )  in  Equation  (2)  shows  that  the  contribution  of  the  pole  to 
the.  field  decreases  exponentially  with  lz+z0l.  In  this  example  lr+r0l  =  2  to  emphasize  the 
near  field  behavior  associated  with  the  pole.  For  larger  values  of  lz+z0l  the  pole  contribution 
would  be  considerable  less.  Equation  (2)  can  be  used  to  estimate  the  magnitude  of  the  pole  con¬ 
tribution  if  the  pole  location,  kr  ,  and  I  z  +Zo!  are  known. 
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CHAPTER  V: 


THE  INVERSION  PROCEDURE 

V.l)  Overview 

I*  Chapter  n  we  developed  the  properties  of  the  Hankel  transform.  These  properties  pro¬ 
vided  our  foundation  for  the  development  of  an  accurate  procedure  to  numerically  generate  syn¬ 
thetic  pressure  fields,  presented  in  Chapter  IV.  In  this  chapter  we  will  use  the  results  of  Chapter 
II  to  explore  the  problem  of  determining  die  plane  wave  reflection  coefficient  from  measure¬ 
ments  of  the  pressure  field  arising  in  a  horizontally  stratified  environment  in  response  to  a  CW 
point  source.  The  estimation  of  the  plane  wave  reflection  coefficient  from  measurements  of  the 
field  is  an  extremely  important  problem.  Determining  the  plane  wave  reflection  coefficient  is  an 
essential  step  in  the  inversion  of  pressure  field  data  to  obtain  the  parameters  of  the  bottom.  In 
this  context  it  is  of  interest  to  geophysicists  and  others  who  wish  to  determine  the  composition  of 
the  ocean  bottom.  The  plane  wave  reflection  coefficient  is  also  used  as  a  geometry  independent 
characterization  of  the  bottom.  As  such,  if  it  is  estimated  in  a  region  from  one  set  of  acoustic 
measurements,  then  the  fields  associated  with  an  arbitrary  source-receiver  geometry  in  that 
region  can  be  determined.  This  is  of  great  value  in  problems  of  acoustic  imaging. 

The  inversion  procedure  that  we  consider  in  this  chapter  was  proposed  by  Frisk, 
Oppenheim,  and  Martinez  [1  ].  It  requires  as  input,  ooherent  measurements  of  the  pressure  field 
as  a  function  of  range  resulting  from  a  CW  point  source  over  a  horizontally  stratified  ocean  bot¬ 
tom.  From  this  the  (complex)  reflected  pressure  field,  Pg  (r ),  is  extracted.  The  Hankel 
transform  of  this  field  is  computed  to  provide  the  depth-dependent  Green’s  function  as  a  func¬ 
tion  of  horizontal  wave  number:1 

a 

G(kr,zj o)  -  fP*(r ]J(dkrf  )rdr  (1) 

0 

Finally,  the  plane  wave  reflection  coefficient  is  obtained  by  multiplying  the  Green's  function  by 


Wswfl 


UspO  dspmdset  Orwa's  function'  to  "Ore* n’»  function'. 


terms  which  compensate  for  the  source  spectrum  and  for  the  source-receiver  distance: 


r(*,)  -  ~i  Vki  - k}e  +IolC  (k,  ,z  ,z0)  (2) 

This  entire  procedure  is  summarized  in  Figure  (V.1.1). 

In  this  chapter  we  will  concentrate  on  the  estimation  of  the  depth-dependent  Green’s  func¬ 
tion.  We  divide  the  issues  addressed  directly  into  the  categories  of  source-field  subtraction,  sam¬ 
pling,  windowing,  and  source-height  variation.  The  issue  of  source-field  subtraction  arises 
because  the  plane  wave  reflection  coefficient  is  directly  related  to  the  reflected  pressure  field  and 
not  the  total  pressure  field,  which  is  measured.  The  issue  of  sampling  covers  the  effects  caused  by 
having  the  pressure  field  available  for  computation  only  on  a  discrete  set  of  points.  We  discuss 
both  the  effect  that  sampling  rate  has  on  the  estimate  of  the  depth-dependent  Green’s  function 
and  the  practical  problem  of  interpolation,  which  is  required  to  move  the  field  from  one  grid  to 
another  (often  to  compensate  for  missing  data  points).  We  develop  a  phase  unwrapping  pro¬ 
cedure  that  allows  us  to  interpolate  the  magnitude  and  unwrapped  phase,  which  vary  slowly 
compared  to  the  quadrature  components. 

In  the  section  on  windowing  we  discuss  the  effect  that  having  the  pressure  field  available 
only  to  a  finite  range  has  on  the  estimate  of  the  depth  dependent  Green’s  function.  We  deter¬ 
mine  the  range  over  which  the  data  must  be  known  in  order  to  accurately  determine  the  depth- 
dependent  Green’s  function.  We  do  this  by  using  the  properties  of  the  Hankel  transform 
developed  in  Section  (II.  6). 

In  the  section  on  source-height  variation  we  exploit  the  results  of  Section  (II. 6)  once  again, 
but  this  time  we  use  them  to  determine  the  effect  that  variations  in  the  source-height  has  on  the 
estimate  for  the  depth-dependent  Green’s  function.  Such  variations  are  inevitable  during  the 
acquisition  of  real  data.  We  illustrate  these  effects  by  considering  the  effect  of  three  specific  vari¬ 


ations. 
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Figure  V.1.1  The  inversion  procedure  to  estimate  the  plane  wave  reflection  coefficient  from  the 
total  field  generated  by  a  CW  point  source  over  a  horizontally  stratified  bottom 


V.2)  Source-field  subtraction 


In  this  section  we  consider  the  removal  of  the  source  field  when  the  source- receiver 
geometry  is  known.  In  Figure  (V.1.1)  we  showed  the  source  field  removed  before  the  Hankel 
transform.  We  did  this  because  conceptually  we  wish  to  deal  with  the  reflected  field  alone.  In  this 
section  we  show  that  numerically  it  is  better  to  remove  the  source  field  in  the  transform  domain, 
after  the  Hankel  transform  of  the  total  field  has  been  computed. 

Because  the  Hankel  transform  is  a  linear  operator,  in  principal  the  estimate  for  the  Green’s 
function  can  be  made  by  subtracting  the  source  field  either  before  transforming: 


G{kr,z,z  o)  =  / 
0 


M')-T 


ttoVrJ+fi-i,,)* 


Vr2+(z  -z0)2 


J^kTr)rdr 


(1) 


or  by  subtracting  in  the  transform  domain: 


G(*r,z,z0)  =  {PT(r)J0(krr)rdr  -  f*  t  ~-J^krr)rdr 

o  oVrz+(z-zay 

which  becomes  upon  performing  the  second  integral  analytically: 


(2) 


G(*r,z,z„)  =  / PT(r)J0(krr)rdr  -  — jA -  (3) 

0  V*o2  -kr2 

If  /’f(r)  is  available  only  over  the  finite  range  0<r<rmix  then  the  field  integrals  can  only  be 
carried  out  to  rmu.  Substituting  rmn  for  in  Equations  (1)  and  (3)  will  make  these  two  formu¬ 
lations  no  longer  equivalent  because  the  analytically  performed  integration  is  not  windowed. 

The  function  transformed  in  Equation  (1)  is  the  reflected  pressure  field,  PK  ( r  ).  In  Section 
(IV. 3)  we  argued  that  the  reflected  field  decayed  asymptotically  as  -^-.1  If  the  total  field,  Pr(r), 

decays  asymptotically  faster  than  — ,  we  can  expect  that  the  formulation  of  Equation  (3)  to 

suffer  leu  from  windowing  effects  than  the  formulation  of  Equation  (1).  We  will  now  show  that 
the  total  field  does  in  fact  decay  faster  than  the  reflected  field  alone.  In  fact,  by  transforming  the 


1)  la  the  ibt«&cs  of  tripped  modes. 
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total  pressure  field  and  then  subtracting  the  source  contribution  in  the  Hankel  domain,  we  have 

performed  the  dual  to  the  removal  of  the  — ...  1 -  singularity  discussed  in  Section  (IV. 3a). 

v  *0 

Before  we  begin,  an  analogy  to  a  similar  procedure  for  determining  the  Fourier  transform 
of  a  function  known  only  over  a  finite  range  but  with  a  large,  known  constant  offset  might  pro¬ 
vide  some  insight.  If  such  a  function  is  (Fourier)  transformed  directly,  the  offset  will  transform 
to  an  impulse  at  the  origin  which  is  smeared  into  the  rest  of  the  transform.  The  smearing  will 
occur  because  the  transform  is  taken  over  only  a  finite  aperture  (windowing).  The  alternative  is 
to  subtract  the  offset  from  the  function,  transform  the  result,  and  add  an  impulse  (with  a 
strength  which  is  determined  analytically  from  the  known  offset)  to  the  origin.  This  second  pro¬ 
cedure  will  give  superior  results  because  the  transform  of  the  offset  is  not  degraded. 

Transforming  the  reflected  field  alone  is  analogous  to  generating  the  Fourier  transform 
directly  from  the  the  field  with  the  known  offset.  In  the  case  of  the  reflected  field,  however, 
instead  of  a  simple  constant  offset,  the  function  has  a  known  asymptotic  behavior.  It  decays  as 

— .  We  are  about  to  show  that  adding  the  source  field  to  the  reflected  pressure  field  is  analogous 
to  subtracting  the  offset  in  the  Fourier  transform  example.  In  the  Hankel  transform  case  we  are 
actually  considering  it  corresponds  ?o  subtracting  a  term  with  the  same  asymptotic  —  behavior. 

The  difference  will  decay  faster  than  — . 


We  begin  by  considering  the  Green’s  function  associated  with  the  total  field  for  Zq>z 
which  is  given  by: 


The 


GT(kr,z,z0)  = 


V*02  -*,2 


r(*> 


/Vtg  -t,2(T0+z)  +  gtVtj 


(4) 


<V*$ -**(»,-») 

- .  i  — —  term  is  the  source  term.  If  we  rewrite  Equation  (4)  in  terms  of  the  reflec- 

Vk$-k? 


don  coefficient  at  kQ,  it  will  be  more  clear  why  adding  this  term  in  the  transform  domain 
corresponds  to  subtracting  the  asymptotic  behavior  in  the  pressure  domain.  We  must  use  the  fact 
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that  for  all  bottoms  r(£g)  =  —1  (except  the  degenerate  case  T(Jtr )  —  1  ).  We  write  Equation 
(4)  as: 

[r(kr)  -  (5) 


Gr(kr,z ,zq)  ~ 


Vkf-k? 

As  we  discussed  in  Section  (IV. 3)  the  asymptotic  behavior  in  the  pressure  domain  is  determined 
by  the  behavior  at  the  singularities  in  the  transform  domain.  [2  ]  At  kr  =  &q  the  phase  term, 


e  21  *r  equals  1  so  that  unlike  the  Green’s  function  associate  with  the  reflected  field 
alone,  the  numerator  of  the  total  Green’s  function  approaches  zero  as  kr  approaches  the 

— 2  ~2  s*n8ularity  at  ^0-  We  wish  now  to  determine  the  contribution  of  the  singularity  at 
V  kfi  —k, 

kr  =  kQ  in  the  total  Green’s  function  in  order  to  show  that  the  "softening"  introduced  by  the 
addition  of  the  source  term  has  made  the  associated  total  field  more  range  limited.  We  can 
bypass  a  great  number  of  issues  by  instead  considering  the  asymptotics  of  the  simplified  Green’s 
function: 


GilmP“(kr?,:Q)  = 


h- 


r(*0)<? 


-2;  V*02-*r2(z) 


-*rJ(z0-O 


(6) 


Vkf-k; 

By  considering  Equation  (6)  we  exclude  those  issues  associated  with  r(jt,.).  Our  examples 
in  the  synthetic  data  generation  section  showed  that  these  terms  do  not  give  rise  to  terms  which 

decay  as  slowly  as  — . 

Equation  (6)  is  the  Green’s  function  for  a  dipole  and  has  the  known  Hankel  transform: 


tt0VrJ+(z-20)J  f*0VrJ+(j+i0)J 

Pd(r)  =  ^  .  ■ 

Vr2+(z  -z0)2  V r2+(z  +z0)2 


(7) 


It  is  well  known  that  this  field  decays  asymptotically  as  — -  and  that  this  asymptotic  behavior 


begins  more  quickly  when  Zq  is  small  than  when  it  is  large. 


Since  the  total  field  will  be  more  range  limited  than  the  reflected  field,  it  is  better  to 
transform  the  total  field  numerically  and  subtract  the  (analytically  determined)  transform  of  the 
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incident  field  than  it  is  to  subtract  the  incident  field  before  transforming. 


V.3)  Sampling 
a)  Overview 

Typically,  data  is  not  available  on  the  grid  required  for  processing.  In  the  experiment  pro¬ 
viding  the  data  for  this  thesis,  for  example,  the  range  values  for  which  the  data  has  been 
obtained  is  determined  by  the  drift  rate  of  the  boat  and  the  source  away  from  the  receivers.  The 
individual  samples  do  not  occur  exactly  where  we  would  like  them  and  while  the  experiment  was 
designed  to  provide  samples  as  close  to  the  Nyquist  rate  as  possible,  typically  the  number  of  sam¬ 
ples  on  averages  is  less  than  we  would  like.  Finally,  there  are  isolated  cases  of  missing  samples,  a 
reality  of  data  collection.  Section  (V.3b)  discusses  the  issues  associated  with  the  average  sam¬ 
pling  rate.  The  issues  associated  with  the  grid  in  general  are  discussed  in  Section  (V.3c). 


b)  Sampling  rate 

In  Chapter  (II.  7)  we  saw  that  when  f(r)  is  sampled  on  approximately  a  linear  grid  and  the 
transform: 


F(p)  =  ff(r)JQ(pr)rdr  (1) 

0 

is  computed  from  these  samples,  then  /  (r )  must  be  sampled  on  a  grid  at  least  as  fine  as  in 

A 

order  to  correctly  perform  the  transform  for  F(p)  negligible  p >A  .  In  this  chapter  we  consider 
the  transform  of  the  pressure  field,  to  obtain  the  depth-dependent  Green’s  function.  This 
transform  has  the  form. 


G(kr,z+z0)  s  -  ]pR{r)Jo{kfr)rdr  (2) 

\kl~kf  0 

G  ( kr,z  +zo)  's  negligible  for  i,. except  possibly  near  the  poles  of  T(kr)  (for  some 


1)  For  some  small  €>0. 
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small  %)  because  when  kr  >fcQ,  it  decays  exponentially.  Consequently  when  there  are  no  poles  in 

T (kr )  for  real  kr ,  then  the  pressure  field  need  only  be  available  on  a  grid  as  fine  as  with 

A 

A  s  to  accurately  determine  C  (kr  ,z  +  Zq)  in  the  region  0<kr<k0.  If  we  wish  to  obtain 
G  (kr  ,z  ~r0)  in  the  region  where  it  is  exponentially  decreasing  (kr>k 0),  however,  we  must  sam¬ 
ple  fast  enough  to  represent  the  signal  in  that  region  as  well. 

If  a  pole  is  present  in  T(ife,)  at  kr  -  kr.  the  Green's  function  will  be  significant  near  krj 
despite  the  exponential  decay.  If  the  presence  of  the  pole  is  ignored  and  the  field  is  transformed 
on  the  grid  ,  then  the  pole  will  be  aliased  into  the  Green’s  function  at  lower  wave  numbers. 

If  there  is  only  one  pole  present  we  can  write  the  Green’s  function,  G  (Jfcr)  (we  suppress  the 
z  variation)  as: 


G(kr)^G(kr)  + 


g-i 

kr~kr2 


for  kr  >k 


(3) 


The  results  of  Chapter  (II. 7)  show  that  if  the  pressure  field  is  transformed  on  the  grid  -j— ,  then, 
the  aliased  Green’s  function  computed  will  approximately  be  given  in  the  region  0 <kr<k  by: 


G(kr)  =  G(kr)  -  V2k/kr-  W(2k-kr)  =  G(kr)  -  y/lk /kr  - 1  ~  (4) 

(2 k~kr)  -kr> 

so  that  the  Greens  function  at  2kQ~kr.  wifi  be  corrupted. 

If  the  amplitude,  a_j,  is  very  small  (which  would  be  the  case  for  large  source-receiveT 
geometries)  and  some  smearing  is  present  due  to  windowing  (the  field  is  not  measured  out  to 
ranges  where  the  trapped  mode  dominates),  we  may  not  see  the  pole’s  effect  and  it  can  safely  be 
ignored. 

In  general  the  possibility  of  trapped  modes  must  be  considered  before  deciding  upon  a  sam¬ 
pling  rate,  particularly  in  geometries  with  small  source-receiver  heights.  For  such  geometries  it  is 

H  IT 

not  always  sufficient  to  sample  at  — — . 
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c)  Sampling  grid 


When  data  is  not  available  on  the  grid  required  for  processing  we  must  first  interpolate. 
Successful  interpolation  is  possible  only  if  the  signal  is  adequately  represented  by  the  original 
samples.  If  we  know  only  that  our  signal  has  a  Hankel  transform  which  is  negligible  beyond 


some  bandwidth,  A,  then  the  signal  is  adequately  represented  by  samples  on  the  grid  -  for 

C 


C  2:A  and  where  rt=  0,  1,  2,  •  •  •  are  the  ordered  zeros  of  Jo(r).  [2  ]  This  is  true  in 
theory.  In  practice,  if  the  the  samples  are  not  originally  spaced  as  required,  it  may  be  impossible 
to  actually  perform  the  interpolation  onto  another  grid.  If  the  samples  are  only  available  on  the 


grid  —  with  C  <A  ,  then  it  is  not  possible  to  interpolate  without  making  additional  assump¬ 


tions. 


We  will  show  that  for  the  class  of  pressure  fields  examined,  an  additional  assumption  seems 
reasonable.  This  assumption  makes  interpolation  possible  even  when  the  sampling  rate  is  slightly 
too  low.  We  will  assume  that  the  magnitude  and  phase  of  our  pressure  fields  2re  smooth  com¬ 
pared  to  their  quadrature  components.  Figure  (V.3c.  1)  shows  the  magnitude  of  the  pressure  field 
associated  with  a  point  source  over  a  pressure  release  bottom.  Between  calculated  points  the 
curve  varies  so  little  and  so  regularly  that  a  plot  of  the  points  appears  to  be  a  smooth  line.  Fig¬ 
ure  (V.3c.2)  shows  the  result  of  first  subsampling  the  points  plotted  in  Figure  (V.3c.l)  (which 


were  available  on  the  grid  — - )  by  a  factor  of  two,  and  then  interpolating  back  onto  the  original 
grid  using  splines.1  The  differences  between  the  two  curves  are  negligiole. 


We  can  compare  this  successful  interpolation  to  the  result  of  subsampling  the  quadrature 
components,  spline  interpolating,  and  computing  the  magnitude.  The  result  of  this  operation  is 
shown  in  Figure  (V.3c.3).  The  apparently  smooth  line  comes  from  the  subsampled  set  of  values 

which  the  splines  was  constrained  to  match  in  the  quadrature  components.  It  actually  consists  of 

1)  Splines  were  used  because  the  original  grid  is  given  by  \n/A  ,  where  are  tbe  ordered  zeros  of 
J q(x  ).  This  grid  is  uneven  and  makes  other  interpolation  schemes  less  desirable. 


every  other  point  of  the  displayed  curve.  The  surrounding  scatter  is  the  magnitude  of  interpo¬ 
lated  points  supplied  by  the  splines.  Clearly,  splines  did  not  successfully  interpolate  the  quadra¬ 
ture  components. 


Figure  (V.3c.4)  shows  the  phase  of  Figure  (V.3c.l)  computed  by 


e 


0) 


where  Pt  is  the  imaginary  component  of  the  field  and  Pr  is  the  real  component.  The  rapid  vari¬ 
ation  in  9  suggests  that  it  is  not  adequately  represented  by  the  grid  upon  which  it  is  presented.  9 
is  not  the  only  representation  of  the  phase  of  the  pressure  field,  however. 


d)  Unwrapping  the  Phase 

The  phase  displayed  in  Figure  (V.3c.4)  is  the  principal  value  of  the  phase,  often  referred  to 
as  the  wrapped  phase.  The  wrapping  comes  about  because  of  the  ambiguity  concerning  which 
phase  should  be  associated  with  the  quadrature  components.  If  9  satisfies: 

Me19  =  P,  +  i  Pi  (1) 

then  so  must  9  +  2ir m  where  m  is  any  integer,  since 

Aie‘<9  +  Ui  m)  =  Mei9e2"i  m  -  Me16  =  Pr  +  iPt  (2) 

Given  just  Pr  and  Pt  there  is  no  way  to  determine  the  correct  value  of  m .  The  arctan  routine 

used  by  Fortran  follows  the  convention  of  choosing  m  such  that 

— * tt  <  §  *  9  +  2irm  <  it  (3) 

The  output  value  §  is  the  principal  value  of  the  phase,  or  the  wrapped  phase. 

If  the  phase  of  the  pressure  field  were  approximately  increasing  at  a  rate  of  kgR  where 
R  ■  Vr2  +  (z  -  z0)2  and  the  field  were  sampled  at  the  Nyquist  rate  of  l/2Jk0,  then  the  phase 
difference  between  samples  would  be  roughly  it  and  the  wrapped  phase  every  sample  or  two 
would  suffer  a  jump  to  a  different  m  in  order  to  satisfy  the  condition  -it  <  9+2'irm  <  it. 
This  would  obscure  any  underlying  regular  behavior  that  we  expect  from  most  physical 
phenomena.  These  frequent  jumps  are  responsible  for  the  rapid  oscillation  apparent  in  the  phase 
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Figure  V.3c.4  The  prirdpal  value  of  the  phase  of  the  field  associated  with  Figure  V.3c.l 


of  Figure  (V.3c.4).  We  wish  to  compensate  for  the  wrapping  that  takes  place  when  the  principal 
value  of  the  phase  is  generated.  To  do  so  we  must  make  use  of  our  knowledge  of  how  the  phase 
is  varying  from  point  to  point. 


We  conjecture  that  the  phase  of  the  total  field  is  dominated  by  a  component  at  the  water 
wave  number  associated  with  the  dominant  specular  path  and  that  the  remaining  portion  of  the 
phase  is  slowly  varying  compared  to  the  sampling  rate.  We  write 

PT{r)  =  M(r  ),'«*'>  (4) 

where  M  (r  )  and  9(r  )  are  real,  and  whte 

e(r)=M+*(r)  (5) 

where  R  *  VV2  -f  (z—  z0)2.  We  will  call  ktfl  the  modeled  portion  of  the  phase  and  t(r)  the 
residual  phase.  We  are  going  to  show  that  as  long  as  the  residual  phase  is  sampled  fast  enough, 
we  can  reconstruct  the  true  phase. 

In  this  notation  the  difference  in  true  phase  value  from  sample  to  sample  can  be  written: 

9(0  "  »('„-!)  =  *o(*»  ~  *„-l)  +  «(*„)  ‘  (6) 

so  that 

9(0  ~  9(r„_1)  ~  ko(R„  -  R.-j)  =  e(R„)  -  e(R„-j)  (7) 

Precisely  stated,  our  requirement  that  the  residual  phase  be  slowly  varying  compared  to  the  sam¬ 
pling-rate  is: 

|c(Rj  -  €(*„-!)!  <TT  for  all  R„  (8) 

To  unwrap  the  phase  we  first  form: 

k> O  ~  “  *0 ~  R»-l) 

=  9(0  -  m„2it  -  -  k0(R,  -  R,-i) 

“  9(0  -  9(r„_1)  -  2 it(m„  -  mw_x)  ^ 

-  «(*„)  -  €(«.-!) 

from  the  measured  data.  We  now  do  the  unwrapping  by  defining  m0  ■  0  and  picking  the 
integers,  m„,  (n  =  1,  2,  •  •  •  )  sequentially  to  satisfy: 
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Mrn)  -  Hrn- 1)  ~  *<>(*«  -  Rn-i)  ~  2 it (rn„  -  m„  _ ,)  =  \t{Rn)  -  t(Rn  -i)  <-rr  (10) 
and  define  the  unwrapped  phase  to  be: 

8(0  s  Hr*)  -  K)2tt  (11) 

Figure  (V.3d.l)  shows  the  result  of  running  this  algorithm  on  the  phase  of  the  synthetic 
data  with  magnitude  shown  in  Figure  (V.3c.l)  and  wrapped  phase  shown  in  Figure  (V.3c.4). 
The  resulting  phase  is  dominated  by  the  linear  term  I(qR  we  defined  in  our  model.  Figure 
(V.3d.2)  shows  the  residual  phase.  The  smooth  and  small  variation  of  the  residual  phase  over 
the  intervals  |r/i-i  •  rn  j  f°r  « ,  is  a  strong  confirmation  of  our  phase  unwrapping  assump¬ 
tion. 

Figures  (V.3d.3a)  and  (V.3d.3b)  present  the  magnitude  and  residual  phase  of  the  fast  bot¬ 
tom  example  of  Section  (IV. 3a).  For  this  example,  too,  the  residual  phase  is  well  behaved. 

Figures  (V.3d.4a)  and  (V.3d.4b)  present  the  magnitude  and  residual  phase  of  the  slow 

* 

speed  layer  example  of  Section  (IV. 3b).  For  this  example,  too,  the  residual  phase  is  well 
behaved.  The  field  in  this  example  was  shown  to  be  dominated  in  the  far  field  by  the  contribu¬ 
tion  due  to  the  pole  beyond  the  water  wave  number.  The  upward  slope  of  the  residual  phase 
apparent  in  Figure  (V.3d.4b)  reflects  the  fact  that  this  pole  is  contributing  the  dominant  com¬ 
ponent  to  the  phase  (in  the  far  field)  which  is  slightly  larger  that  the  IcqR  term  subtracted  out. 


e)  Interpolating  the  magnitude  and  unwrapped  phase 

In  Figures  (V.3c.l)  and  (V.3c.2)  we  showed  that  the  magnitude  of  the  dipole  field  could  be 

K  K 

up-sampled  from  the  grid  ——  to  .  In  Section  (V.3d)  we  saw  that  the  unwrapped  and  resi- 
dual  phase  components  enjoy  smooth,  regular  variation  ideally  suited  for  interpolation.  Figure 

K  K 

(V.3e.l)  shows  the  residual  phase  for  the  dipole  field  up-sampled  from  the  y  to  the  -y  grid. 


It  is  indistinguishable  from  the  residual  phase  generated  on  the  -y  grid  shown  in  Figure 


(V.3d,2). 
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Figure  V.3d.l  Unwrapped  phase  corresponding  to  Figure  V.3c.4 


MAGNITUDE 


We  now  show  that  for  the  dipole  field  we  can  actually  interpolate  the  magnitude  and 
unwrapped  phase  to  increase  the  sampling  rate  of  the  quadrature  components  of  the  field.  This 
allows  us  to  determine  the  Green's  function  to  a  higher  horizontal  wave  number  than  the 
Nyquist  criteria  applied  to  our  original  sampling  scheme  would  have  us  believe. 


SI- 

We  recall  from  Section  (II.7)  that  if  the  pressure  field  for  the  dipole  on  the  grid  —  were 

transformed  and  displayed  in  the  range  0<fcr<2  the  result  would  be  severely  aliased  and  com¬ 
pletely  inaccurate  in  the  region  Kkr  <2.  To  obtain  a  transform  accurate  on  0<ir<2,  the  qua- 

K 

drature  components  must  be  at  least  sampled  on  the  grid  — We  can  still  obtain  the  transform 

in  the  range  0<Jkr<2,  never-the-less,  by  interpolating  the  field  onto  the  grid  through  its 

magnitude  and  unwrapped  phase.  Figures  (V.3e.2a)  and  (V.3e.2b)  show  the  magnitude  and 
phase  of  the  transform  generated  by  such  a  procedure.  First  the  magnitude  and  residual  phase  of 

p  ( —)  were  generated.  These  were  up-sampled  as  just  discussed.  From  this  up-sampled  magni¬ 
tude  and  residual  phase  (and  the  modeled,  k(fl  ,  portion  of  the  phase)  the  associated  quadrature 
components  were  generated.  This  was  transformed.  Figures  (V.3e.3a)  and  (V.3e.3b)  show  the 

,  ,  K  s 

magnitude  and  phase  of  the  Hankel  transform  of  p(-y)  generated  without  interpolation.  Only 

small  differences  in  the  magnitude  are  apparent.  The  phase  curves  also  display  only  small  differ¬ 
ences  though  in  the  inhomogeneous  region  (where  the  phase  is  oscillating  rapidly  as  evidenced  by 
the  two  parallel  lines)  the  small  difference  has  caused  a  slightly  different  picture  of  the  oscilla¬ 
tions.  By  contrast.  Figures  (V.3e.4a)  and  (V.3e.4b)  present  the  magnitude  and  phase  of  the 


Hankel  transform  of  p(y-)  up-sampled  by  direct  spline  interpolation  of  its  quadrature  com¬ 


ponents.  Clearly,  once  again,  a  direct  interpolation  of  the  quadrature  components  did  not  work. 

We  apply  this  scheme  for  interpolating  the  magnitude  and  residual  phase  to  the  field  of  the 
fast  bottom  example  of  Section  (IV. 3a).  We  first  generate  the  magnitude  and  unwrapped  phase 


Figure  V-3e.4b  Phase  of  Hankel  transform  of 
been  up-sampled  through  its  quadrature 
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.  \n  x„ 

of  p  interpolate  up  to  the  grid  generate  the  quadrature  components,  and  then 

J  Zi 

Hankel  transform.  The  magnitude  and  phase  of  the  result  is  shown  in  Figures  (V.3e.5a)  and 

(V.3e.5b).  The  magnitude  and  phase  of  the  correct  transform  (of  p(~~)  generated  without 

using  this  interpolation  scheme)  is  shown  in  Figures  (V.3e.6a)  and  (V.3e.6b).  We  see  that  the 
Hankel  transform  of  the  up-sampled  data  and  the  Hankel  transform  of  the  data  originally  avail¬ 
able  on  the  fine  grid  do  not  agree  exactly.  Figures  (V.3e.7a)  and  (V.3e.7b)  present  the  magni¬ 
tude  and  phase  of  their  complex  difference  and  Figure  (V.3e.8)  presents  the  magnitude  of  the 
Hankel  transform  of  that  complex  difference.  This  transform  represents  the  errors  made  in  the 
pressure  domain  by  our  up-sampling  procedure  that  gave  rise  to  the  error  in  the  Green’s  func¬ 
tion.  We  see  that  practically  all  the  error  energy  was  concentrated  at  the  origin.  This  error  could 
be  due  to  a  breakdown  in  our  phase  unwrapping  assumption  near  the  origin  or  to  a  poor  han¬ 
dling  of  the  rapid  change  in  magnitude  by  the  splines.  This  problem  can  be  corrected  by  a  dense 
sampling  of  the  original  field  near  the  origin  so  that  there  is  no  room  for  interpolation  error 


f)  Phase  unwrapping  errors 


At  this  point  we  consider  briefly  the  kinds  of  error  that  might  be  expected  when  the 
assumption  underlying  this  phase  unwrapping  technique  is  violated.  If  for  some  n 


I »(*.)  -  -  *o (*„  “  >  *  (1) 

the  wrong  m„  will  be  chosen.  From  that  point  on,  each  mt  (Jt  =  n ,  n  + 1,  •  •  •  )  chosen  by  the 
procedure  will  also  be  wrong  by  the  same  amount  A  plot  of  this  error  is  a  step  function  of 
height  fh„  -  m,  centered  at  n  as  shown  in  Figure  (V.3f.l).  If  multiple  violations  occur,  the 
error  will  look  like  the  sum  of  step  functions  as  illustrated  in  Figure  (V.3f.2).  The  smoothness 
apparent  in  the  residual  phase  in  all  of  our  examples  suggests  that  no  errors  have  occurred. 

If  the  phase  unwrapping  scheme  is  used  to  interpolate  the  field,  these  errors  are  not  serious. 
As  part  of  the  interpolation  procedure,  the  quadrature  components  are  regenerated  from  the 
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Figure  V.3e.5a  Log-magnitude  of  Hankel  transform  of  the  field  for  the  fast  bottom  example 
after  it  has  been  up-sampled  through  its  magnitude  and  residual  phase 


Figure  -V.3e.6a  Log-magnitude  of  Hankel  transform  of  the  held  for  the  fast  bottom  example 
generated  originally  on  a  fine  grid 
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Figure  V.3e.6b  Phase  of  Hank  el  transform  of  the  field  for  the  fast  bottom  example  generated 
originally  oa  a  fine  grid 


PHASE 


interpolated  magnitude  and  phase.  The  error  curve  shown  in  Figure  (V.3f.2)  would  have  no 
affect  on  the  quadrature  components  regenerated  from  the  unwrapped  phase.  In  general,  after 
the  interpolation,  the  error  curve  will  not  be  a  simple  sum  of  steps  but  will  be  smeared  by  the 
interpolator.  This  will  usually  affect  the  quadrature  components.  If  the  interpolator  is  well 
chosen,  the  leakage  will  be  small  and  limited  to  the  area  near  the  error.  Finally  we  note  that 
errors  in  the  phase  unwrapping  scheme  will  occur  when  the  unmodeled  portion  of  the  phase  is 
varying  rapidly  between  samples.  When  this  happens  the  interpolator  is  likely  to  have  difficulties 
even  without  errors  in  the  unwrapped  phase  and  this  scheme  is  probably  not  appropriate. 


V.4)  Windowing 

In  Section  (II. 6)  we  stated  that  in  terms  of  resolution  the  Hankel  transform  behaves  very 
much  like  a  Fourier  transform.  We  wish  to  consider  the  resolution  required  to  generate  the 
Green’s  function  and  the  window  that  this  implies. 

The  total  Green’s  function  is  given  by 


Gr(*r,z)  =  -  ,  J.  -  -  [r(*,)g'V**~tf2(j+ro)  +  <r'Vi<?~*,2(l',o)]  when  z>z0  (1) 

V  k§  —kf 

The  most  rapid  variations  in  Gj  (excluding  possible  poles  in  the  reflection  coefficient  beyond  the 


water  wave  number)  occur  near  kr  =  ko ■  When  r(kr  )  is  smooth  compared  to  — 


Vk#-k2 


^Vki-kfr+io) 

the  rapidity  of  these  variations  is  dominated  by  the  - ,  ■ —  term.  With  a  windowed 

V*02-*r2 

sample  of  the  pressure  field  we  can  not  hope  to  determine  the  exact  behavior  of  G  at  kr  -  k  0, 

when  z  +zq  is  large  the  rapid  variation  in  G  is  due  to  primarily  the  e*  k^~k'^r  +I°^  term.  We 
can  obtain  an  ad  hoc  estimate  of  the  resolution  we  require  by  considering  the  lobe  widths  associ¬ 
ated  with  the  phase  for  kr  near  k0. 


That  is  we  define  K*  by  the  relation: 


*  fifj  ’tHMP . 


-  166  - 


'l 


f, 


and  define  the  lobe,  8,  by: 


Vkf^k^(z+z0)  * 

8*  “  ^r,i i-l  ~^r,n 


nit 


(2) 


0) 


If  we  use  (Jko“^r^i)(^0+*rji)  =  (-^—  )2  then  when  kr  .  ~  *  we  have: 

z  +z0 


and 


8.  * 


(*o 

Ka)  *  2*0 

ntt 

Z 

z+zQ 

1 

IT 

(n2-(n-l)2) 

1 

ir 

2*o 

z+z0 

. 

2*o 

z+z  o 

(4) 


(2n  +1) 


Section  (II.6c)  indicates  that  the  required  window  width,  B ,  is  related  to  the  desired  resolution 
approximately  as: 


aa33i  =  3(*+zq)2  1  1 

8  8„  ip’  2*0  2 n  1 


(6) 


Thus  to  resolve  the  lobe  doses  to  k  when  z  +zq  =  136  and  kQ  -  .9246  we  require  a  win¬ 
dow  of  about: 

\2 


B 


=  _ 1 _ 

\  it  J  (2)(.9246) 


=  3*  (1013)  =  3040  meters 


(7) 


V.5)  Source-Height  Variation 
a)  General  expression 

The  procedure  proposed  to  estimate  the  plane  wave  reflection  coeffident,  r(ir),  and 
shown  in  Figure  (V.1.1),  requires  that  the  pressure  field  be  measured  with  the  source  at  a  fixed 
height,  zo.  [1  ]  Frequently,  experimental  conditions  cause  the  source-height  to  vary.  In  this  sec¬ 
tion  we  will  explore  the  effect  that  a  varying  source-height  has  on  the  estimate,  f  (*f  ). 

Instead  of  considering  the  effect  of  a  varying  source-height  on  the  estimate  for  the  plane 
wave  reflection  coeffident  directly,  we  will  consider  its  effect  on  the  depth-dependent  Green’s 
function  given  by:1 

1)  W#  will  suppress  the  dependence  of  Hi*  depth-dependent  Orsen't  function  oa  z  sad  Zq. 


mm 


>*•***-:*•  5^  *-• 


r 
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_/L  \  fn  /  /.  \  J  1  )  i  Vt^  -*,2(i  +!o) 

G(*r)  =  J>j?(>-^VQ(fc,r)r^r  =  —7 - . je 

0  V  kn  —kr 


' r/  jo  “v  —  '  v*o2-*; 

which  is  the  Hankel  transform  of  the  reflected  pressure  field.  By  considering  the  effect  on  the 
depth-dependent  Green’s  function  we  can  make  use  of  the  properties  of  the  Hankel  transform 
that  we  derived  in  Chapter  II.  The  plane  wave  reflection  coefficient  is  determined  by  multiplying 
the  estimate  of  the  depth-dependent  Green's  function  by  terms  which  compensate  for  the  source 
strength  and  the  source-receiver  separation  as  was  shown  in  Equation  (V.1.2). 

We  consider  the  effect  of  a  source  height  given  by 

2o(r)  =  2o+*(r)  (2) 

To  explore  the  effect  of  Hankel  transforming  a  pressure  field  measured  at  a  source  height  that  is 

a  function  of  range,  we  write  the  Green’s  function  estimated  by  Hankel  transforming  this  field 


» 

{r  a  {r))Jol*,r)rdr 

-i|i» 


■foUrttd?  Uo(*rr)rdr 


We  now  define: 


»(*,>$)  -  if elV^hir)J 0(ir )J 0(kr r )rdr 


which  with  (3)  becomes: 


=  /G(OH(*r,$)dS  (5) 

0 

Equations  (4)  and  (5)  exactly  describe  the  effect  that  source-height  variation  has  on  the 
estimate  of  the  depth-dependent  Green’s  function.  As  written,  however,  they  do  not  provide 
much  insight  into  what  variations  are  tolerable  or  into  the  qualitative  effect  of  source-height  vari¬ 
ation.  To  provide  us  with  this  insight  we  develop  an  approximate  expression  for  H  (kr  ,§)  by 


! 
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using  the  windowing  result  of  Section  (11.6).  We  do  this  by  considering  e  as  a  win 

dow.  The  result  will  be  reasonable  provided  that  the  Fourier  transform  (in  r)  ol  e  is 

narrow,  as  discussed  in  Section  (II. 6b). 


We  write 


H(kr,Q  3  ^Jei,“h(r)J 0(ir)J Q(krr)rdr  with  o>  =  V*<?-42 


Hkr~0 

The  Hankel  transform  of  Jq(&)  equals - - - so  that : 


vr,H(kr,o  =  i -  VW*-  -«> 


where 


W&,)  -  J 


This  provides  us  with  approximate  expressions  for  the  kernel,  H(kr,£),  and  an  approximate 
expression  for  the  estimated  Green’s  function  in  terms  of  the  actual  Green  s  function. 

<?(*,)  - 

In  the  following  sections  we  apply  this  result  to  some  special  cases. 


b)  Particular  variations 


i)  No  source-height  variation 


When  the  source-height  is  constant,  li(r)  s  0.  For  this  case  our  approximate  result  above 
gives  W£kr -i)  =  2ir8(*r  -4)  and  G  (*,)  =  G  (*r),  which  is  as  we  would  expect. 

11)  Linear  source-height  variation 

If  the  source-height  varies  linearly  then  h(r)  =  ar  and 
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with  a)  s  a  y/icQ—f?. 


To  perform  the  integral  in  Equation  (2)  we  expand  the  exponential  in  its  Taylor  series  to 
obtain: 


w5(*r)  =  J  2 


H . 


e  '  dr 


~x/i  =0  n  * 


n =0  n  •  -* 


(3) 


=  2it£-iI^S(*r-na) 

*=0  B ! 

We  use  this  to  determine  the  effect  of  the  real  cosinusoidal  variation: 


h(r)  = 


ae,ar  +  ae~lar 


(4) 


Substituting  (4)  into  Equation  (V.5a.8),  for  cosinusoidal  variation  W ^(kr)  is  seen  to  be  given  by: 

W,/*,)  =  Je  2  e  2  e'Vdr  (5) 

-X 

Equation  (5)  is  the  Fourier  transform  of  a  product  of  terms  in  the  form  of  Equation  (2).  Conse¬ 
quently  we  can  write  Equation  (5)  as  the  convolution  of  terms  in  the  form  of  Equation  (3): 


n  “0  "! 

IcVkf^ 


*«(»,)  - 

Zir  *„0  n ! 


kr 


2rr  2  h(kr+na) 


ii  “0 


n! 


X  X 

=  2tt2  2 

a  -Obi  =0 


(6) 


n\m ! 


5(ikr  -(n  -m)a) 


If  we  perform  the  integration  (V.5a.9)  we  obtain: 


i  i  f 

■  ’hi  0  I 

kr>{n-m)S~ 


taVtj  -(*,-(n  -m)q)2  \ 

-  - V*r-(n  -m)aC  (kr - (n  —m )a) 


•  Ui! 


(7) 


The  cosinusoidal  source-height  variation  with  an  amplitude,  a ,  and  a  frequency  a,  has  the 
effect  of  reverberating,  or  comb  filtering,  \^k~rG{kr)  in  two  dimensions.  The  impulses  of  the 


required  (or  this  case. 
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filter  are  spaced  at  the  variation  frequency  a.  The  .  weightings,  which  we  define  to  be 
w  (kr  ]tn  ,n  ),  determine  the  envelope  of  the  reverberation.  They  depend  on  the  amplitude,  a : 


w  (ir  ;m  ,n  )  = 


2 


m+n  iiaV^T 


F 


We  can  write  the  estimate  for  the  Green's  function  in  terms  of  these  weighting  functions  as; 
»  » 

x  22 

G(kr)  ~  kr>(n-m)  a  W(K  ~(n  -m)a,m,n)Vkr  ~(n -m  )a  G  (kr  ~(n  -m  )a)  (9) 

The  weighting  functions,  w  (kr  -  (n  —m  )cr/rt  ,n),  are  greatest  when  m  -  n  =  nmtx  and 
decay  rapidly  from  that  point  in  m  and  n .  This  result  can  be  shown  by  replacing  the  factorials  in 

(8)  with  Stirling’s  approximation  (excellent  even  for  small  n :  n !  ~  ~n  j  )  and  defining 
ia  -k~ 

x  =  - - - .  The  weighting  functions  then  become: 


(?f 


w(kr.m  ,n)  =  ~  —  =  — yl -  —  —yl — (—  1  (10) 

ml  n\  ^2-nm  m  J  ^2tui  (  n  j 

term  has  its  maximum  at  m  =jr-  and  falls  off  in  m  with  greater  than  geometric 


decay.  The  —y  term  pulls  this  maximum  only  very  slightly  lower. 

’  m 

The  result  is  that  w  (kr  /n  ,n  )  is  large  for  m  ,  n  Vito  — k,2  and  small  elsewhere.  When 
|  Vito2  -k}a  s  ■  1  «1  we  can  ignore  the  (n-m  )a  term  in  (9)  so  that  d  (Jtr)  is  given 


approximately  by 


By  defining: 


,  il  m 


v7j|  kr>{*-m)* 


'k,-(n  -m)aC  (*,-(«  -m)a) 
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C(*r*)-(i)"  X 

m  =max(Q,-n) 


[-ir 


we  can  obtain  the  result,  valid  for  j-^—  «  1  that 


(n  -rm )! 


Vk^G(kr)  =  2  C(kr*)Vkr-n<xG(kr-na)  (13) 

»  -  -X 

A  perspective  plot  of  C(kr,n )/(/)."  is  presented  in  Figure  (V.Sc.la)  for  the  case  a  =  3 
and  Jt0  =  .9246159.  The  back  of  this  figure  corresponds  to  kr  =  k0  and  consequently  kz  =  0. 
The  Green’s  function  in  this  region  corresponds  to  plane  wave  components  of  the  field  which  are 
directed  entirely  in  the  radial  direction  and  which  do  not  vary  in  z .  Figure  (V.Sc.lb)  presents 
the  slice  of  Figure  (V.5c.la)  corresponding  to  this  region,  C  (k0,n).  C  {k0,n)  is  zero  everywhere 
except  at  n  =0,  where  it  is  1.  Referring  to  Equation  (13)  we  see  that  the  degraded  estimate  of 
the  Green’s  function  at  kr  =  to  is  given  by: 


G(*0)  =  2  C(k0,n)Vk0-naG(k,-na)  (14) 

Substituting  for  C  (io»n)  Equation  (12)  we  see  that 

G(*o)  =  C(*0)  (15) 

The  portion  of  the  spectrum,  kr  =ko,  corresponds  to  field  components  that  do  not  vary  in 
z .  It  is  reasonable,  then,  that  the  cosinusoidal  source-height  variation  did  not  affect  that  portion 
of  the  angular  spectrum . 

In  Figure  (V.Sc.la)  moving  forward  towards  the  leading  edge  corresponds  to  decreasing  kr 
and  increasing  kt .  With  decreasing  kr,  C(kr,n)  becomes  increasingly  less  impulsive,  indicating 
greater  amounts  of  degradation.  Figure  (V.Sc.lc)  presents  the  slice  C(0,n).  This  slice 
corresponds  to  that  portion  of  the  angular  spectrum  which  has  the  maximum  amount  of  vertical 
variation.  In  Figure  (V.Sc.lc)  the  value  C(0,0)  is  not  even  as  large  as  the  adjacent  values, 
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Figure  V.Sc.lc  C(0,/i)  with  a  =  3 
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C  (0,-1)  and  C  (0,1)-  The  Greens  function  will  be  degraded  by  cosinusoidal  source-height  vari¬ 
ation  in  this  region. 

Figure  (V.5c.2)  presents  a  perspective  plot  of  C (kr ,n)/(i)n  for  the  case  a  =  12  and 
kQ  =  .9246159.  Once  again  C (ko,n)  is  the  discrete  delta  function,  S(n),  and  the  Green’s  func¬ 
tion  will  not  be  degraded  at  kr  =  jfc0.  Because  a  is  larger  now,  \  k%-~k}  of  Equation  (12) 

grows  more  rapidly  as  kr  becomes  smaller  than  it  did  for  a  =  3.  Asa  result  the  figure  shows 
that  serious  degradation  begins  for  k,  much  closer  to  ko-  The  increased  amplitude,  a ,  has 
resulted  in  an  increased  amount  of  degradation.  The  product,  ak,  =  aW k <j  ~k?,  determines 
the  severity  of  this  effect. 

We  note  also  that  because  of  the  (i)*  factor  in  Equation  (12),  the  phase  of  C(kr,n) 


ically  affect  the  phase  of  the  estimated  Green’s  function,  (j  (kr),  even  before  it  significantly 
affects  the  magnitude. 

Thus  we  have  seen  that  the  effect  of  sinusoidal  source-height  variation  is  to  comb-filter  the 
estimate  of  y/k^G(kr).  The  spacing  between  impulses  in  the  comb  filter  is  the  frequency  of  the 
source-height  variation.  The  amplitude  of  the  source-height  variation  and  the  vertical  wave 
number,  k£  —  fcr2,  together  determine  the  weightings  of  the  impulses.  When  the  product  of  the 
amplitude  and  the  vertical  wave  number  is  small,  the  only  contribution  comes  from  the  low  lag 
components.  As  this  product  increases,  the  higher  lag  components  begin  to  contribute  and  the 
comb-  filtering  will  become  increasingly  apparent.  If  the  frequency  of  the  source-height  variation 
is  very  low,  causing  the  spacing  of  the  impulses  in  the  filter  to  be  very  small,  the  degradation 
will  appear  as  a  smearing. 

V.S)  Summary 

In  this  chapter  we  have  studied  the  issues  associated  with  the  inversion  of  pressure  field 
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data  through  the  Hankel  transform  to  estimate  the  depth-dependent  Green’s  function  and  the 
plane  wave  reflection  coefficient.  We  have  developed  a  phase  unwrapping  procedure  that  allows 
us  to  interpolate  the  magnitude  and  unwrapped  phase  and  thereby  determine  from  the  set  of 
field  samples  available,  the  values  of  the  field  at  the  ranges  we  require  for  processing.  We  have 
also  shown  that  it  is  better  to  estimate  the  total  depth-dependent  Green’s  function  from  the 
Hankel  transform  of  the  total  field,  and  to  later  remove  the  affects  of  the  source.  Finally,  we 
have  examined  the  effects  of  source-height  variation  to  help  us  understand  the  possible  degrada¬ 
tion  that  this  effect  would  would  introduce  into  the  depth-dependent  Green’s  function  estimated 
from  real  data. 

We  are  ready  to  perform  a  preliminary  processing  of  real  data. 
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CHAPTER  VI) 

Inverting  Real  and  Realistic  Data 

In  Chapter  V  we  described  the  procedure  for  inverting  coherent  field  measurements  arising 
in  response  to  a  CW  point  source,  to  obtain  the  plane  wave  reflection  coefficient.  In  that 
chapter  we  addressed  some  of  the  practical  issues  that  must  be  faced  when  real  data  is  to  be 
inverted.  In  this  chapter  we  perform  a  preliminary  inversion  of  real  data.  [1  ]  To  help  interpret 
the  results,  in  parallel  we  invert  data  generated  synthetically  for  a  realistic  geometry  and  set  of 
bottom  parameters. 

The  real  data  that  we  invert  was  obtained  by  G.  Frisk,  J.  Doutt,  and  E.  Hays  in  1981. 
The  associated  experimental  geometry  was  described  in  Section  (1.6)  and  is  presented  again  in 
Figure  (VI. 1).  We  will  be  using  the  data  obtained  from  the  lower  receiver  shown  in  this  figure. 
In  Figure  (VI.2)  we  present  a  velocity  profile  and  density  parameters  for  a  bottom  that  we 
believe  is  comparable  to  the  bottom  where  the  real  data  was  taken.  We  use  this  geometry,  velo¬ 
city  profile,  and  these  density  parameters  to  generate  the  synthetic  data  of  this  chapter.  This 
synthetic  data  is  generated  using  the  hybrid  procedure  described  in  Chapter  IV  and  the  numeri¬ 
cal  Hankel  transform  that  was  described  in  Section  (ID. 7)  -1  [2  ]  The  efficiency  of  this  Hankel 
transform  algorithm  made  it  possible  to  obtain  high  quality  results  over  a  large  range  that  would 
otherwise  not  have  been  practical. 

We  begin  by  generating  the  synthetic  data  for  this  geometry  and  bottom.  We  use  the 
numerical  procedure  described  in  Section  (IV.  2)  to  generate  the  plane  wave  reflection  coeffi¬ 
cient,  r(*r),  as  a  function  of  horizontal  wave  number.  Its  magnitude  and  phase  are  presented  in 
Figures  (VI.3a)  and  (VI.3b).  We  see  that  a  pole  is  present  in  r(&r)  beyond  the  water  wave 
number.  This  pole  is  due  to  the  low  speed  channel  just  below  the  water-bottom  interface. 
Because  the  source -‘-receiver  height  is  large,  this  pole  will  contribute  an  insignificant  amount  to 

1)  This  algorithm  was  implemented  in  Fortran  on  a  VAX-1 1/780  by  Mike  Wengrovitz. 
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Figure  VI.l  Experimental  geometry  associated  with  the  real  data 
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Figure  VI. 2  The  velocity  profile  and  density  parameters  for  a  bottom  comparable  to  that  where 
the  real  data  was  taken 
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Figure  VL3b  Phase  of  the  reflection  coefficient  calculated  for  bottom  parameters  shown  in  Fig* 
ure  VI.2 
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the  pressure  field  and  need  not  be  removed  as  v/as  done  in  Section  (17. 3b).  Consequently  we 
can  generate  the  field  with  the  hybrid  procedure  described  in  Section  (IV. 3a).  The  magnitude 
and  residual  phase  of  the  associated  synthetic  field  are  presented  in  Figures  (VI.4a)  and  (VI.4b). 
In  these  figures  very  little  high  frequency  ripple  is  apparent  even  at  large  ranges,  implying  that 
the  field  is  indeed  adequately  represented  and  not  suffering  from  spatial  aliasing. 

Figures  (Vl.Sa)  and  (Vl.Sb)  present  the  magnitude  and  residual  phase  of  the  synthetic 
field  after  inclusion  of  the  incident  field.  The  regular  behavior  in  these  plots  suggests  that  the 
magnitude  and  residual  phase  are  good  representations  of  the  total  field.  As  further  confirma¬ 
tion  of  the  validity  of  the  total  synthetic  fields  generated  for  this  example,  we  present  the  output 
of  a  ray  program  that  was  run  for  this  profile  in  Figure  (VI. 5c). 1,2  The  two  synthetic  fields  axe 
in  good  agreement  except  in  the  region  of  the  caustic,  1500m  <  r  <  2000m ,  where  the  ray 
method  is  known  to  be  inaccurate. 

Figures  (VI.6a)  and  (VI.6b)  present  the  magnitude  and  residual  phase  of  the  real  data 
(which  includes  the  source  field),  fit  the  region  beyond  the  first  hundred  meters,  the  magnitude 
and  residual  phase  of  the  real  data  behave  regularly,  which  gives  us  confidence  in  them.  The 
interference  pattern  apparent  in  the  magnitude  is  similar  to  that  of  the  synthetic  data.  The  zeros 
in  the  magnitude  are  well  matched  by  the  the  changes  in  the  residual  phase  for  large  ranges.  The 
first  few  hundred  meters  of  the  residual  phase,  however,  looks  significantly  different  from  the 
residual  phase  of  the  synthetic  field.  In  this  region,  changes  in  the  source-height  have  their 
greatest  effect  on  the  measured  field  because  the  geometry  is  most  significantly  affected  by 
source-height  variation  in  this  region.  We  recall  that  the  residual  phase  is  given  by: 

e(r)  =  8(r)  -  t0Vr2+(r  -z0)2  (1) 

The  large  negative  slope  of  the  residual  phase  for  low  ranges  could  be  due  either  to  an  estimate 

of  J k<)  which  was  too  large,  in  which  case  the  residual  phase  would  display  a  negative  phase 

everywhere,  or  to  an  estimate  of  (z  -z0)2  which  was  too  large.  We  believe  that  this 

1)  Bat  with  •  ihghtty  different  loans  height  of  125  meters  rather  than  135  meters. 

2)  1  with  to  dunk  Jim  Doutt  tor  providing  this  synthetic  field. 
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Figure  VI.5e  Log-magnitude  of  the  total  field  calculated  by  a  ray  method 
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uncharacteristic  behavior  is  due  to  imperfect  knowledge  of  the  source  height  in  that  region.  The 
gentle  negative  slope  of  the  unwrapped  phase  for  the  real  data  for  large  ranges  is  probably  due 
to  a  slight  overestimate  of  the  water  wave-number,  leg. 


Before  we  attempt  to  invert  the  real  experimental  data  to  estimate  the  depth-dependent 
Green's  function  two  major  factors  must  be  considered.  First,  the  experimental  data  is  available 
only  over  a  finite  range  and  second,  it  is  available  only  at  discrete  points  which  are  not  spaced 
properly  for  our  processing.  The  first  issue  can  be  resolved  by  referring  to  Chapter  V  where  we 
showed  that  for  the  source-height  and  geometry  used  to  obtain  the  experimental  data,  it  was  only 
necessary  to  know  the  field  out  to  about  3040  meters  to  minimize  the  degradation  due  to  win¬ 
dowing.  The  experimental  data  is  available  to  6000  meters.  We  believe,  therefor,  that  window¬ 
ing  should  not  prevent  its  successful  inversion.  The  second  issue  can  also  be  resolved  by  refer¬ 
ence  to  Chapter  V  where  we  showed  that  by  interpolating  the  magnitude  and  unwrapped  phase 
it  was  often  possible  to  translate  the  pressure  field  data  available  on  one  set  of  ranges  to  another. 
We  will  use  the  procedure  developed  there  to  interpolate  the  experimental  data  onto  the  set  of 
ranges  that  we  require  for  processing  by  the  Hankel  transform.  In  parallel  we  will  process  the 
synthetic  data.  The  processed  synthetic  data  provides  a  useful  measure  of  the  success  of  our  pro¬ 
cessing  because  the  depth-dependent  Green’s  function  that  we  obtain  can  be  compared  to  the 
true  depth-dependent  Green’s  function  which  is  known  for  the  synthetic  data,  and  presented  in 
Figures  (VI. 7a)  and  (VI.7b). 

Figures  (VI.8a)  and  (VI.8b)  present  the  magnitude  and  phase  of  the  Green’s  function  cal¬ 
culated  by  processing  the  synthetic  data.  The  synthetic  data  was  originally  available  on  the  grid 
air  it  =  0,  1,  2,  ■  •  ■  .  It  was  linearly  interpolated  (through  its  magnitude  and  unwrapped 


phase  as  described  in  Section  (V.3e))  onto  the  grid  required  for  processing,  —  for 

A 


n=0,l,2,  •••  with/4  =  1.2  and  where  n=0,  1,  2,  •••  are  the  zeros  of  /o(x)- 


The  agreement  between  the  estimate  of  the  synthetic  Green’s  function  obtained  by  process¬ 
ing  the  synthetic  field  and  shown  in  Figures  (VI.8a)  and  (VI.8b)  and  the  true  Green’s  function 
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for  the  synthetic  data  shown  in  Figures  (VI. 7a)  and  (VI. 7b)  is  particularly  the  agree¬ 

ment  in  the  magnitudes.  The  phases  differ  slightly  for  low  horizontal  wave  numbers  We  believe 
that  this  is  due  to  small  errors  in  the  synthetic  field  for  low  ranges.  The  phases  differ  dramati¬ 
cally  in  the  evanescent  region  beyond  the  water  wave  number,  where  the  magnitude  of  the 
Green’s  function  is  very  small  and  consequently  the  phase  is  probably  dominated  by  noise.  The 
agreement  in  general  between  the  true  Green’s  function  for  the  synthetic  data  and  the  Green’s 
function  estimated  from  the  synthetic  data  is  excellent,  however,  and  confirms  the  results  of 
Chapter  V  which  indicated  that  for  the  sampling  rate  and  range  of  values  over  which  the  data  is 
known,  it  should  be  posable  to  determine  the  depth-dependent  Green’s  function. 


Figures  (VI.9a)  and  (VI.9b)  present  the  magnitude  and  phase  of  the  Green’s  function  cal¬ 
culated  from  the  real  data.  Except  for  low  wave  numbers,  the  magnitude  of  this  Green's  function 
displays  many  of  the  features  of  the  synthetic  Green’s  function,  including  the  same  overall 


envelope  due  to  the 


L  2  _  L  2 

kr 


source  spectrum  term,  and  the  interference  pattern  arising  from 


the  interaction  of  that  portion  of  the  Green’s  function  associated  with  the  source  and  that  portion 
associated  with  the  reflected  field.  The  total  Green’s  function  also  decays  rapidly  at  the  water 


wave  number,  as  it  should  due  the  the  e 
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migration  terms.  In  the  evanescent 


region,  kr>  fc0,  we  see  only  noise,  comparable  to  the  noise  we  see  superimposed  upon  the  rest 


of  the  spectrum. 


At  low  horizontal  wave  numbers  the  Green’s  function  for  the  real  data  does  not  look  like 


the  Green’s  function  for  the  synthetic  data.  Very  near  the  origin  we  see  a  large  peak  not 
apparent  in  the  total  Green’s  function  for  the  synthetic  data.  This  peak  is  probably  due  to  con¬ 
centration  of  noise  power  there  by  the  Hankel  transform  as  discussed  in  Section  (II.8).  For 
slightly  larger  wave  numbers  the  magnitude  displays  a  jagged  appearance  not  seen  in  the  total 
Green’s  function  for  the  synthetic  data.  In  this  region,  tne  stationary  phase  approximation  for 
the  Sommerfeld  integral  is  fairly  good,  allowing  us  to  associate  the  behavior  of  the  Green’s  func¬ 
tion  at  low  horizontal  wave  numbers  with  the  behavior  of  the  pressure  field  at  low  ranges.  The 
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Figure  VI. 9a  Phase  of  the  depth-dependent  Green’s  function  estimated  by  taking  the  Hankei 
transform  of  the  interpolated  real  data 
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uncharacteristic  behavior  of  the  Green’s  function  at  low  horizontal  wave  numbers  is  consistent 
with  the  uncharacteristic  behavior  of  the  residual  phase  that  we  observed  for  low  ranges  and  may 
be  due  to  variations  in  the  source-height.  Some  of  this  apparent  jitter  in  the  Green’s  function  of 
the  real  data  may  be  due  in  part  to  variation  in  the  source-height.  A  rough  sampling  of  the 
source-height  over  the  course  of  the  experiment  was  available  from  the  experimental  records. 
We  interpolated  between  available  samples  using  splines  to  obtain  a  rough  estimate  of  the 
source-height  variation  present  during  the  course  of  the  experiment.  The  result  is  presented  in 
Figure  (VL 10).  This  curve  is  sufficiently  similar  to  the  sum  of  the  two  low  frequency  cosines  dis¬ 
cussed  in  Section  (V.5b.iii)  to  qualitatively  interpret  the  effect  of  source-height  variation  for  this 
experiment  in  term  of  the  results  presented  there.  The  analysis  of  Section  (V.5)  shows  that 
sinusoidal  variation  in  the  source-height  causes  the  estimated  Green’s  function  to  be  a  rever¬ 
berant  version  of  the  true  Green’s  function,  particularly  for  low  kr  corresponding  to  larpe  kx . 
Because  the  frequency  of  the  variation  is  very  small,  the  main  effect  is  to  smear  the  estimate  of 
VTrG  (kr).  As  stated  in  that  section,  the  phase  of  the  estimated  Green’s  function  might  be  more 
seriously  corrupted  than  its  magnitude.  The  phase  of  the  depth-dependent  Green’s  function 
estimated  from  the  real  data  and  shown  in  Figure  (VI.9b)  does  not  strongly  resemble  the  phase 
of  the  synthetic  Green’s  function.  The  overall  good  appearance  of  the  magnitude  of  the  total 
Green’s  function  and  the  poor  appearance  of  it  phase  is  consistent  with  the  the  degradation  that 
would  be  expected  from  source-height  variation. 

Figures  (VI.  11a)  and  (Vl.llb)  show  the  magnitude  and  phase  of  the  plane  wave  reflection 
coefficient  generated  from  the  Green’s  function  calculated  from  the  synthetic  data  and  shown  in 
Figures  (VI.8a)  and  (VI.8b).  Figures  (VI. 12a)  and  (VI.12b)  present  the  magnitude  and  phase 
of  the  plane  wave  reflection  coefficient  calculated  from  the  Green's  function  for  the  real  data. 
The  estimate  for  the  reflection  coefficient  for  the  real  data  does  not  appear  to  be  a  good  one  at 
this  time. 

Because  the  plane  wave  reflection  coefficient  is  obtained  from  the  total  depth-dependent 
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Figure  VI. 11a  Magnitude  of  the  plane  wave  reflection  coefficient  generated  from  the  depth- 
dependent  Green’s  function  estimated  from  the  synthetic  data 
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Figure  VI. 12a  Magnitude  of  the  plane  wave  reflection  coefficient  generated  from  the  depth* 
dependent  Green’s  function  estimated  from  the  real  data 
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Green’s  function  by  first  coherently  subtracting  the  source  contribution  and  then  multiplying  by  a 

term  with  a  rapidly  varying  phase  (<?  1  0  '  ‘  2°'  )  errors  in  the  phase  of  the  total  Green  s 

function  would  seriously  degrade  the  estimate  of  plane  wave  reflection  coefficient.  The  estimate 
for  the  reflection  coefficient  is  probably  much  worse  than  the  estimate  for  the  Green’s  function 
because  of  the  phase  errors  in  the  estimate  for  the  total  depth-dependent  Green’s  function. 


In  conclusion,  we  believe  most  of  the  error  apparent  in  the  Green’s  function  for  the  real 
data  to  be  due  to  variation  in  the  source- height  near  the  origin.  Direct  evidence  of  this  is  the 
anomalous  residual  phase  variation  in  the  region  r  <  300  meters.  The  error  in  the  estimated 
Green’s  function  for  very  small  horizontal  wave  numbers  is  probably  due  to  addidve  noise.  The 
errors  in  the  estimated  reflection  coefficient  are  probably  due  to  imperfect  knowledge  of  the 
source-receiver  geometry  that  affects  the  coherent  additions.  Overall,  however,  we  are  greatly 
encouraged  by  the  good  appearance  of  the  magnitude  of  the  total  depth-dependent  Green’s  func¬ 
tion  determined  from  the  real  data.  The  interference  structure  and  the  overall  envelope  suggest 
that  we  are  very  dose  to  being  able  to  estimate  the  plane  wave  reflection  coefficient  from  real 
data.  Work  still  needs  to  be  done  to  compensate  for  the  effect  of  source-height  variation. 


Hie  potential  returns  from  the  successful  inversion  of  pressure  field  data  to  obtain  the  plane 
wave  reflection  coefficient  are  enormous.  Such  a  successful  inversion  is  a  vital  step  in  the  process 
of  inferring  the  physical  parameters  of  the  bottom  from  acoustic  measurements.  [3, 4  ]  The  abil¬ 
ity  to  make  such  inferences  is  of  great  interest  to  oceanographers  and  to  exploration  geophysi¬ 
cists.  A  successful  inversion  would  also  make  it  possible  to  predict  the  fields  associated  with  an 
arbitrary  source-receiver  geometry  from  one  set  of  measurements.  This  would  greatly  facilitate 
acoustic  imaging  in  the  ocean. 
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CHAPTER  VII: 

CONTRIBUTIONS  AND  FUTURE  WORK 

Vn.l)  Contributions 

In  this  thesis  we  have  studied  both  the  numerical  generation  of  synthetic  pressure  fields 
from  the  plane  wave  reflection  coefficient  and  the  inversion  of  measured  pressure  field  data  to 
estimate  the  plane  wave  reflection  coefficient.  We  developed  and  implemented  algorithms  that 
efficiently  generate  high  quality  synthetic  fields.  We  studied  the  major  issues  affecting  the  inver¬ 
sion  of  experimental  data  and  were  able  to  estimate  the  depth-dependent  Green's  function  from 
measured  data  taken  in  the  ocean  with  a  high  degree  of  success.  We  isolated  source-height  varia¬ 
tion  as  a  major  factor  preventing  the  successful  estimation  of  the  plane  wave  reflection  coefficient 
at  this  time. 

As  a  foundation  for  our  studies  we  explored  the  Hankel  transform  in  depth.  In  Chapter  II 
we  derived  a  number  of  important  properties  including  the  effects  that  windowing  and  sampling 
a  function  have  on  its  Hankel  transform.  Our  sampling  results  show  that  the  associated  degrada¬ 
tion  is  often  a  more  severe  problem  for  the  Hankel  transform  than  for  the  Fourier  transform,  hi 

particular  it  can  seriously  degrade  synthetically  generated  pressure  fields  which  decay  as  —  or 

r 

even  more  slowly  and  its  effect  should  always  be  carefully  considered. 

In  Chapter  II  we  also  studied  the  noise  properties  of  the  Hankel  transform.  We  showed 
that  if  a  function  is  sampled  on  a  square  root  grid  in  a  noisy  environment,  its  Hankel  transform 
will  have  superior  noise  properties  more  characteristic  of  the  underlying  two  dimension  Fourier 
transform  which  the  Hankel  transform  represents  in  the  presence  of  cylindrical  symmetry. 

In  Chapter  in  we  considered  a  number  of  numerical  techniques  for  performing  the  Hankel 
transform.  We  presented  new  results  strengthening  existing  procedures  such  as  the  asymptotic 
and  back  smear  methods  as  well  as  an  efficient,  exact  method  developed  as  part  of  this  thesis. 

In  our  development  of  algorithms  to  generate  high  quality  synthetic  data  we  presented  a 
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number  of  hybrid  numerical-analytical  techniques  that  greatly  improve  the  quality  of  synthetic 
data.  In  the  course  of  developing  a  technique  that  can  adequately  handle  the  effects  of  guided 
modes  in  slow  speed  layers  under  the  ocean  bottom,  we  derived  an  expression  that  may  be  of  use 
for  developing  held  expressions  in  modal  expansions  required  to  be  accurate  in  both  the  near 
and  far  fields.  We  also  presented  a  well  behaved  numerical  procedure  for  implementing  the 
Thomson-Haskell  approach  for  generating  the  plane  wave  reflection  coefficient. 

In  Chapter  V  we  developed  the  major  issues  affecting  the  inversion  of  measured  field  to 
obtain  the  plane  wave  reflection  coefficient.  On  the  basis  of  this  development  we  were  able  to 
identify  the  sources  of  error  in  an  actual  inversion.  The  phase  unwrapping  and  interpolation 
results  presented  in  this  chapter  also  significantly  improved  the  results  of  the  processing  of  the 
experimental  data  in  Chapter  VI. 

In  Chapter  VI  we  performed  a  preliminary  inversion  of  real  data  to  obtain  estimates  for 
the  depth-dependent  Green’s  function  and  the  plane  wave  reflection  coefficient.  The  results 
presented  in  this  chapter  represent  a  significant  advance  towards  the  complete  inversion  of  meas¬ 
ured  pressure  field  data  to  obtain  the  plane  wave  reflection  coefficient.  We  were  able  to  generate 
a  good  estimate  for  the  depth-dependent  Green’s  function  and  were  able  to  associate  the  effects 
of  source-height  variation  with  the  degradation  in  the  estimate  for  the  plane  wave  reflection 
coefficient 

At  this  point  work  is  continuing  towards  the  complete  estimation  of  the  plane  wave  reflec¬ 
tion  coefficient  from  real  data.  The  foundations  laid  by  the  work  presented  in  this  thesis  provide 
a  strong  base  for  future  work  in  this  area.  In  addition  they  suggest  research  in  a  number  of 
related  areas.  Some  of  these  are  presented  in  the  next  section. 


VII. 2)  Future  work 
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a)  Cylindrical  to  Cartesian  Coordinate  Systems 

In  this  thesis  we  have  dealt  with  problems  cast  in  a  cylindrical  coordinate  system.  In  that 
coordinate  system  the  familiar  Fourier  transform  of  cartesian  systems  became  the  less  familiar 
Hankei  transform.  In  that  coordinate  system  the  dean  properties  of  additive  white  Gaussian 
noise  through  the  Fourier  transform  were  obscured  obscured  until  a  square  root  grid  was  intro¬ 
duced.  In  that  coordinate  system  windowing  and  aliasing  approximately  affected  ^rf  (r) 


instead  of  /  (r).  In  that  coordinate  system  the  familiar  impulse  8(jr)  became 


Mil 


In  that 


coordinate  system  the  — r  operator  became  V2  »  — —  + - so  that  die  operator  which 

fix2  Hr2  r  ar 


nulls  a  pole  in  a  cartesian  coordinate  system: 


a 2 


+  0? 


ax*  '  p' 

became  the  less  familiar: 


/(*)  = 


ax 3 


Pi 


/  -j~eivxdp  =  /  — j — 'j-e^dp  =  -8(x) 
_  nr— nr  _  nr— n* 


p2+p2 


-* p  -pf 


-*  p  -pf 


fv2  +  p 2 1/ (r )  =  fv2  +  P/2]/  2l  2J<tpr)?dp  =  J— 3 — ~~Jo(f>r)pdp  =  (2) 

l  J  L  Jop2-p,  0  Pz-P/  r 

In  short  we  frequently  found  that  familiar  problems  in  a  cartesian  system  became  more  difficult 
when  caste  in  a  cylindrical  coordinate  system.  The  reverse  is  also  true,  however.  In  Section 
(III.  7)  we  developed  an  efficient  numerical  algorithm  for  the  Hankei  transform  by  mapping  it 
into  a  Fourier  transform.  The  mapping  was  accomplished  with  the  Abel  transform: 


(1) 


(3) 


The  Abel  transform  also  serves  to  map  other  linear  operators  in  cylindrical  form  into  linear 
operators  in  cartesian  form  (it  must  do  so  for  any  function  that  can  be  represented  by  a  Hankei 

,  d2 

transform).  It  maps  V*  into  — -  for  example,  in  the  sense  that: 


8x 


‘M-£4w) 


(4) 


i 


In  Section  (III. 7)  we  developed  an  efficient  numerical  algorithm  for  evaluating  the  Abel 
transform.  One  exciting  area  of  future  research  is  the  extension  of  maximum  entropy  and  other 
spectral  estimation  techniques  into  the  cylindrical  domain  as  we  now  describe. 

An  estimation  scheme  that  might  be  of  value  for  estimating  the  plane  wave  reflection  coef¬ 
ficient  is  illustrated  in  Figure  (VII. 2a.  1).  Instead  of  estimating  the  plane  wave  reflection  coeffi¬ 
cient  directly  we  estimate  the  position  and  residue  of  its  poles  (and  possibly  zeros)  in  the  complex 
plane.  We  do  this  because  the  estimation  of  parameters  instead  of  a  function  is  a  much  better 
posed  problem  when  the  signal  available  for  analysis  has  been  corrupted.  Because  the  plane  wave 
reflection  coefficient  is  related  to  the  measured  field  by  the  (cylindrical)  Hankel  transform  and 
not  the  Fourier  transform,  spectral  estimation  techniques  available  in  the  literature  of  digital  sig¬ 
nal  processing  do  not  apply.  If  we  first  process  the  pressure  field  data  with  an  Abel  transform, 
however,  the  resulting  signal  has  the  same  poles  and  zeros  but  now  in  its  Fourier  transform. 
Modern  spectral  estimation  techniques  can  therefore  be  used.  We  introduce  the  caution  that  the 
effect  that  branch  cuts  have  on  this  procedure  must  be  studied  with  care. 

b)  Analytical-Numerical  Algorithms 

The  hybrid  analytical-numerical  technique  used  to  implement  the  Abel  transform  in  Section 
(III.6)  is  a  very  general  procedure  and  springs  from  classical  numerical  methods  of  long  standing. 
Traditionally,  difficult  integrals  are  evaluated  numerically  by  removing  their  singular  behavior  as 
much  as  possible  through  coordinate  changes  and  changes  of  variables  and  then  numerically 
transforming  the  result.  The  success  of  the  hybrid  method  points  out  that  in  fact  it  is  often  desir¬ 
able  to  do  just  the  opposite.  The  integral  should  be  manipulated  to  produce  as  much  singular 
behavior  as  possible.  The  singularities  can  be  integrated  analytically  and  will  not  suffer  from 
numerical  degradation.  If  the  singularities  are  removed  properly,  the  remaining  numerical  por¬ 
tion  of  the  integral  will  be  well  behaved  where  it  dominates  and  subordinate  to  the  analytically 
determined  portions  of  the  integral  where  it  does  not.  The  art  in  this  procedure  is  casting  the 
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analytically  determined  portions  of  the  integral  carefully  to  insure  that  the  numerical  portion 
does  not  have  infinities  or  undesirable  asymptotics  to  cancel.  Had  the  e  ~br  factor  not  been 
included  in  the  procedure  of  Section  (HI. 6),  for  example,  the  hybrid  method  for  the  Abel 
transform  would  not  have  worked. 

Manipulating  functions  so  that  they  can  best  be  represented  by  parameterized  functions  and 
samples  is  part  of  the  general  issue  of  computer  representation.  As  software  systems  become 
smarter  this  kind  of  approach  will  become  increasingly  more  important. 


c)  Waveguides 

In  the  course  of  generating  synthetic  data,  we  evaluated  the  integral: 


/(r,z,p,  )  s  j*  ,7*'-" r  T'l*'  *2  p2U^o(Pr)prfp 
0  V*2_p2  p2-p2 


and  showed  that  it  satisfied 


(1) 


[Vr2+Pi2j/(r,z,Pi)  -  f  r-  -  =  o  (2) 

We  associated  f(r,z,P/)  with  the  contribution  of  the  pote  at  Pi  because  our  integrals 


always  included  the  terms  — >  ■  —  c1  **  as  well.  The  advantage  of  this  formulation 

V*2_p2 


was 


that  7(r  ,z  ,p()  is  everywhere  finite,  even  at  r  =0.  The  classical  contribution  associated  with  a 
pole  is 


when  Im(P<)>0  (3) 

The  Hankel  function  above  has  a  logarithmic  singularity  as  the  origin.  Physically,  poles  in 
the  depth-dependent  Green’s  function  make  only  finite  contributions  to  the  field.  The  migration 
term  *<V^5-p2^ '  windows  the  pole  in  the  Hankel  domain  so  that  its  contribution  to  the  pressure 
field  is  everywhere  finite.  For  this  reason  we  included  the  migration  term  into  our  pole  expres¬ 
sion.  For  convenience  we  also  included  the  source  term  — ;  ■  -  ■  . 

V*2—  p2 
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Our  formulation  was  carefully  constructed  to  insure  that  the  numerics  were  not  required  to 
generate  infinities.  It  is  potentially  useful  for  many  problems  other  than  those  considered 
directly  in  this  thesis.  One  set  of  problems  concerns  the  calculation  of  fields  inside  (possibly 
leaky)  waveguides. 

We  develop  an  expression  for  the  field  inside  a  (dielectric)  waveguide  arising  from  a  point 
source  using  the  plane  wave  formulation  of  this  thesis.1  Figure  (VII.2c.l)  presents  the  geometry 
of  the  waveguide  and  the  waves  present. 

The  radial  and  time  variation  of  all  fields  is  given  by: 


JQ(krr)e -f“»  (4) 

and  will  be  suppressed.  We  will  use  fi  for  the  vertical  wave  number.  (3  and  kr  axe  related 

through 


k2  =  tf+k2  (5) 

The  source  field  is  given  by  Ptel^x  ro1  and  is  the  portion  of  the  field  that  would  be 
present  even  in  the  absence  of  impedance  contrasts.  We  specify  the  boundary  conditions  at 
z  =  ft  by  giving  the  plane  wave  reflection  coefficient  there,  TT(kr),  and  at  z  =  0  by  T B(kr). 
These  two  interfaces  together  give  rise  to  an  up-going  wave,  P  +e‘^z ,  and  a  down-going  wave 
P  ~e~'^2 ,  that  would  not  be  present  without  the  impedance  contrasts  in  the  regions  zSrft  and 
zk 0.  The  total  up-going  and  down-going  fields  for  Zq<z  t£ft  is  given  by: 


+  Pjemz~^  UP 
P-e’1*2  DOWN 

-  The  plane  wave  reflection  coefficient  at  z  =ft  provides  the  boundary  condition: 

P-e-'?h 

T  "  Py*  +  P/e'p(*",°) 

In  the  region  O^z  <Zq  we  have 

P-e’'*2  +  DOWN 

P+e1*2  UP 


(6) 

(7) 

(8) 


A  derivation  for  the  field  inside  such  e  waveguide  can  also  be  found  in  [1  ] 
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n  i  pO  a 

r  =  _ +e _ _  (9) 

B  P-e~‘W  +  Pjt f,pi°  P_  +  p,e^z° 

Pj  is  the  known  point  source  strength  -^-^====-  =  FT  and  TB  are  given  as  well. 
The  two  boundary  conditions  are  sufficient  to  determine  P  +  and  P  _  from  these  quantities.  We 


IYP+e'P*  +  rTP,ei?he  ,pr°  =  P_e~'P* 

rBp _  +  r Bp,e‘^  =  p  + 

Solving  for  P  +  and  P  -  we  have 

r -  r7Vp(*~Io)l/>, 

P+  =  A - 7 - s  C+P, 

-^"'P*  -  rre'P* 


= 


r^ciP(*+l0)  _  eim-i o)  p 

-J-e-’P*  -  rBe'Pfc 

rr 


The  total  field  in  the  waveguide  (  0<z  <h  )  is  given  by 

*('■*)-{  vir 7?  I*  +c-e 


.V*2-*A 


]^0(*rO*r 


(13) 


c+  = 


-i-tf-‘P*  -  rre'P* 

1  A 


r  #P(^+Jo)  _ 

1  J»C  c 

C-m-~ -  (15) 

-J-e-'P*  -  rBe*'P* 
lr 

The  zeros  of  e-,Ph  -  T B  r^e'  P*  =  0  contribute  poles  to  the  depth-dependent  Green’s 
function  and  give  rise  to  the  modes  of  the  waveguide.  Each  of  these  poles  in  Equation  (13) 
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makes  a  contribution  to  the  field  of  the  form 

X  _ 

I(r,z,pt)  =  J0(pr)pdp  (16) 

o  V*2_p2  P  “Pi 

and  the  development  we  used  for  calculating  fields  in  the  presence  of  poles  applies. 

It  should  be  noted  that  the  inverse  problem,  that  of  resolving  modes  in  a  waveguide,  can  be 
cast  into  the  classical  signal  processing  form  of  finding  poles  in  the  Fourier  transform  by  first 
generating  the  Abel  transform  of  the  pressure  field!  The  effect  of  branch  lines  on  this  approach 
needs  to  be  studied,  however. 

It  is  also  possible  to  construct  nulling  operators  to  estimate  the  pole  positions  (as  is  done  in 
Maximum  Entropy  spectral  estimation  techniques)  by  using  Equation  (2),  which  does  indude 
some  branch  line  effects. 
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APPENDIX  1: 

DETERMINING  LIMITS  FOR  THE  ALIASING  RESULT  OF  SECTION  (H.7) 
In  this  appendix  we  determine  the  limits  as  N-&  for  the  expressions: 


sin  N  itx  itx 

- cos  — 

itx  4 

sin  — 

2 


cos  N  rrx  itx 
sm  2 

To  evaluate  the  limit  in  Equation  (1)  we  use  the  known  limit: 


sm2irx  (N  -  — )  . 

..  2  sin  2irxN  v 

lim - ; - -  lim - cos  irx  =  >8(x  — * )  (3) 

v-.*  smirx  /v-*  sm  -nx  * 

The  limit  in  Equation  (3)  is  zero  except  at  the  zeros  of  smirx ,  which  occur  at  x  -  k  for 
k  =  0,  ±1,  ±2,  •  •  •  .  At  these  points  coswx  =  cosit*  =  (-1)*.  The  effect  of  deleting  the  cosirx 
term  is  to  generate  alternating  signs: 


lim  ”n2Tr/V*  =  Y(-1)*S (x-*) 
smirx  j 


A  change  of  variables  shows  that 

X  .  L 

sm^y 

We  determine  the  limit  of  Equation  (1)  by  multiplying  by  00a— •  to  find: 

sinirNx  ...  itx  vx-^,  ,^„,x  ,  ^ 

“  T  T  ~  6(j~k) 

sin  it * 

2 

We  can  simplify  Equation  (6)  by  using  the  fact  that  /  (x  )8(x  -x0)  =  /  (x0)8(x  -x0)  and  that 

1  it  =0,4,8,  •  •  • 

IT*  0  *=1,5,9,  •  •• 

“*  2  ~  -1*  =2,6,10,-  •• 

0  *=3,7,11,  •  •  • 


to  find: 
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lim 


*a7,N*co*3j-  =  2(-l)*8(y  -  2k) 

sirnr^-  * 

2 


This  is  the  result  needed  in  the  text. 


(8) 


To  evaluate  the  limit  of  Equation  (2): 


cos  N  irx  irx 

- sm  — 

irx  4 

an  -r— 


(9) 


we  first  consider  the  behavior  of 


cosNirx 


near  x  =  0.  The  behavior  at  the  other  zeros  of  sin-?- 


sin- 


TTX 


will  be  similar.  We  assume  that  x  is  sufficiently  small  so  that  sin-^-  can  be  replaced  the  first 


term  of  its  Taylor  series,  We  consider: 


cosA'irx 
lim  - 

W-*  TTX 
2 


(10) 


Instead  of  evaluating  the  limit  directly,  we  look  at  the  Fourier  transform  of  the  limit.  We 


evaluate  this  by  taking  the  limit  of  the  Fourier  transforms  of  each  term  and  write:1 


lim 

cos/V  irx 

—  lim  FT 

'  ' 

COsWlTX 

TTX 

iV-x 

TTX 

i 

2 

* 

2 

(11) 


The  Fourier  transform  of  C0S^1TX  can  be  found  by  convolving  the  Fourier  transform  of 

T 

2 

cosAf  irx  with  the  Fourier  transform  of  — . 

irx 


so  that: 


FT 


cosN-nx 


irx 

2 


-"{"""HI-} 


(12) 


l)Tbete  step*  an  be  rigorously  justified  by  using  generalized  functions.  [1  ]  They  presume  that  two  func¬ 
tions  are  said  to  be  equal  if  the  result  of  convolving  their  difference  with  any  bud-limited  function  is  al¬ 
ways  zero.  Alternately,  two  functions  are  said  to  be  equal  if  the  Fourier  transform  of  their  difference  is 
zero  for  any  finite  band. 


FT 


cosAf  TTX 

ITX 

2 
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»(/“)  +5(/+f)'* 


[-2*  sgn  (/  )  ] 


This  function  is  plotted  in  Figure  Al.l.  It  is  given  by 


FT 


oosN  rrx 


TTX 

2 


A'  N 

4l  *  ~~2 

0 

— 4*  f  <f 


(13) 


(14) 


As  N-x i  this  Fourier  transform  becomes  equal  to  zero  over  any  finite  interval.  Consequently,  as 
a  generalized  function 


i  • 


V 


Hm  cosNjrx_  ^  q 
nx  * 
2 


(15) 


where  0(  means  0  as  a  generalized  function.  Basically  this  means  that  as  N  -=°  CW^~'1WC  oscillates 

ITX 

2 

rapidly  around  zero  in  such  an  manner  that  when  it  is  convolved  with  any  bandlimited  function, 
the  result  is  zero.  Since  we  will  only  use  this  functions  inside  integrals  (strictly  speaking  impulses 
are  only  defined  inside  integrals)  we  will  simply  call  it  zero. 

t  COS N  *1TX 

Given  that  this  limit  is  zero,  the  limit  of  - ,  which  is  simply  a  periodically  repealed 

ITX 

smT 

version  of  cosNl,x  near  the  zeros  of  sin(-^~ )  (with  sign  changes),  must  also  be  zero.  Since 
itx  2 


sin”-  is  finite  everywhere,  we  have  the  result  needed  in  the  text  by  multiplying: 


..  cosAfirx  .  itx  n 

Hm - sin—  =  0 

.  itx  2 

nn__ 


(16) 


t 
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APPENDIX  II: 

THE  VALUE  OF  THE  KERNEL 

FOR  THE  NUMERICAL  PORTION  OF  THE  HYBRID  ALGORITHM 
AT  THE  WATER  WAVE  NUMBER 

In  this  appendix  we  derive  the  value  of  L(k q)  discussed  in  Section  (IV. 3a).  £.(*,)  is 
defined  by  Equation  (IV. 3a.  15)  as: 


,  1  L(k')  L(k°)  iVki-k,1  III 


and  we  seek  to  evaluate  the  limit 


lim  L  (kr) 

t,-*o 


under  the  condition  that  the  impedance  of  the  bottom,  Z  (*r),  is  finite  at  kr  —  ij. 

At  kr  =  *0  Equation  (1)  takes  the  indeterminate  form  We  evaluate  the  limit  (2)  by 
using  L’  Hospital’s  rule: 

i{r(jkr)  -  r(fc0)V -^-«fr(*,)  -  r(k0)V 

ii _ ' _ l _ -  is™  — l_i - 1 -  ™ 


'kl  ~kr2 


—  lim 


After  separating  out  the  terms  that  approach  zero  as  k,-*o  this  expression  becomes: 


.  -iVfco-kr2  „  iy/ki-k} |. 

hm  L(kr)  =  lint  - - - i(kr)e 

A_-H in  Km 


We  now  express  r(*,)  in  terms  of  the  characteristic  impedance  of  the  upper  half  space,  Zo, 
and  the  impedance  at  the  interface  which  we  will  denote  at  Zj.  Both  Zo  and  Z i  are  functions  of 


kr  in  general.  In  terms  of  these  r(i,)  is  given  by: 


r(*,)  = 


Zi~Z0 

Zj+Zo 


Taking  derivatives  we  find: 


*>  - 2 


-  224  - 


We  now  use  Che  characteristic  impedance  of  the  upper  half  space: 


-k,2 


where  po  is  the  density  of  the  half  space,  o>  is  the  temporal  frequency,  and  k0  *  where  Co  is 

the  speed  of  sound  in  the  upper  half  space.  Substituting  Equation  (7)  into  Equation  (6)  and 
evaluating  20  we  have: 


PQU  trP|)U  _  *rPo“  _ 

po*>  pfc3  zm -* 3>+2 2lP0^Vki -k}  +  p w 

z'  i re 

Substituting  (8)  into  (4)  we  find 


_  *rPO“  „ 

|W’  *  v^?z' 


.  —  ViK- 

*,-*«  (  r  *,-*0  K  Z?(*0  -*r2)+2Z1p0o)Vfc02  -V+p0V 


L(*o)  = 


-2iZi(*o) 


provided  that  lim  "S /kj—  kr2Zt(kr)  =  0 

*r~*0 


If  the  interface  is  between  two  isovelocity  half  spaces,  the  expression  (10)  for  L  (*o)  can  be 

Pi<i> 

written  directly  in  terms  of  the  material  parameters.  For  this  case  Zj(*r)  =  -f'\  2-  zi(*o)  «* 

\kx-kr 

finite  because  fe1  #to  ( if  Jtj  =  ko  there  would  be  no  interface).  L  (to)  is  given  by 


*■(*<>)  = 


Po  y/kf-kH 


(11) 
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APPENDIX  III: 

EVALUATION  OF  THE  POLE  CONTRIBUTION  TO  THE  FIELD  FOR  SECTION  (IV.3B) 
Here  we  evaluate  the  pole  contribution  to  Equation  (1): 


/ (r ;pf )  =  S~T~2T7r  e/V*1~*’i|z  l/0(pr )p<i p 

0  PZ~P f  V*2_p2 


(1) 


We  evaluate  Equation  (1)  by  determining  a  partial  differential  equation  that  it  satisfies  and  solv¬ 
ing  that  equation. 

Taking  the  second  partial  derivative  of  Equation  (1)  with  respect  to  z  we  have 

a£/(r,z;Pi)  = 

If  we  use 


8(z)/0O  =  8(z)/(0)  for  any  /(z) 


(3) 


then  Equation  (2)  becomes 


£rnr,z-Pl)  =  / p2 — p2  VJv g<Vrr?1, l^pr)prfp-26(z)/-^/0(pr)pdp  (4) 


Putting  it  all  together  we  have: 


(Pi2  “  *2)  /(r,z;p()  =  /— rj — -e'V*2-^l,l/0(Pr)pdp  -  25 (z)  /  2  1  2  Jo^pdffi) 
riz  J  0  Vl!-p f  o  P  “  P) 

If  we  define  p2  ss  p2  —  k2  ,  choosing  real  part  of  P  >  0,  and  use 


*  i  jVTTI,  i 

Jt7,2  ,g  ~p  lllJ0(pr)pdp  =  — 7- 
0  Vjfc2_  2  Vr 


2+z2 


(6) 


together  with 


/”2~2"/°^pr)prfp  =  Ko(~s8"Um(Pi)]  tPi')  =  -=^/fi1)(p,r)  when  >  0  (7) 

0  P  P(  i 


then  (A2.5)  becomes 


il-B2 
2  P 


flz 


I  git  Vf*+i2 

/  (r  ,r  ;p^ )  =  — - — —  +  i  (p,r) 

r^+z* 


(8) 


when  Im  (P/ )  >  0. 


* 
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The  Green’s  function  for  this  differential  equation  is  given  by 

G(r’z*5)  =  (9) 
Using  this  for  the  impulsive  response  and  convolving  with  the  continuous  driving  function  to 
obtain  the  particular  solution  we  obtain: 

=  -TJ-/  -7"'''  -  ^■ffi1)(p,r)e_pl*1  (10) 

Ip  *  Vr2+.e2  ■ip 


When  Im  (p, )  >  0. 


The  general  form  of  this  expression  which  is  valid  for  all  p/  is  given  by: 


I(r,z;pi)=  2p~/  %r-2"e  WpiO*  31,1  (11) 
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